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Abstract  of  Dissertation  Presented  to  the  Graduate  School 
of  the  University  of  Florida  in  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of  Doctor  of  Philosophy 

ASSOCIATION  MODELS  FOR  CROSS-CLASSIFICATIONS 
HAVING  ORDERED  CATEGORIES 

By 

Abbas  Kezouh 
August  1984 

Chairman:  Alan  G.  Agresti 

Major  Department:  Department  of  Statistics 

This  dissertation  considers  the  analysis  of  categorical  data  that 
includes  ordinal  variables.  Three  specific  problems  are  examined: 

(1)  generalizing  the  models  introduced  by  Goodman 

(2)  testing  for  independence  between  ordered  row  categorical 
variables  and  column  categorical  variables  of  an  rxc  cross- 
classification 

(3)  estimating,  in  two-way  tables,  association  parameters  under 
the  constraint  that  they  are  monotonic. 

Goodman's  models  are  generalized  so  that  they  can  be  applied  to 
cross-classifications  having  arbitrary  dimensions.  Primary  emphasis  is 
devoted  to  three-way  cross-classifications.  This  includes  modeling, 
interpretation,  and  fitting  of  the  models. 

Using  a simpler  model  than  the  one  proposed  by  Haberman,  we  pro- 
pose test  statistics  for  independence  in  an  rxc  cross-classification. 


vi 


After  deriving  some  properties  possessed  by  the  proposed  statistics, 
we  present  the  results  of  a simulation  study  to  investigate  the 
performance  of  the  statistics,  comparing  them  to  the  conventional 
Pearson  and  likelihood-ratio  chi-square  statistics. 

Finally,  we  propose  methods  that  produce  monotonic  estimated 
parameters  for  two  models  that  will  be  described  in  the  second  chapter. 
Also,  a simulation  study  will  be  used  to  examine  the  performance  of 
the  estimates  produced  by  the  proposed  method  compared  to  the  estimates 
produced  by  standard  methods. 
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CHAPTER  ONE 
INTRODUCTION 

Extensive  consideration  has  been  given  to  methods  for  analysis  of 
categorical  data  in  the  past  two  decades,  as  evidenced  by  the  volume 
of  work  on  the  subject  published  during  this  period.  In  particular, 
a comprehensive  treatment  has  been  given  by  Bishop  et  al . (1975). 

However,  the  early  texts  seem  to  lack  sufficient  attention  to  ordinal 
categorical  data.  To  overcome  this  shortcoming,  standard  log! inear 
models  were  merely  applied  to  ordinal  settings.  It  is  only  in  recent 
years  that  specialized  models  for  ordinal  data  have  been  investigated. 
This  effort  was  pioneered  by  Haberman  (1974a)  and  Goodman  (1979). 

In  the  standard  logl inear  models,  where  the  variables  are  treated 
as  if  they  were  nominal  in  scale,  the  problem  is  twofold.  On  one  hand, 
nontrivial  models  do  not  exist  to  give  a parsimonious  description  of 
the  table,  in  the  sense  that  the  models  are  complex.  For  instance,  for 
two-way  tables,  in  the  scheme  of  the  nominal  hierarchical  models,  the 
only  model  more  complex  than  the  independence  model  is  the  saturated 
model.  On  the  other  hand,  they  lack  power  for  detecting  certain  associa- 
tions, especially  if  a natural  order  applies  to  at  least  one  of  the 
categorical  variables.  If  we  take  into  account  the  apparent  ordering 
of  some  or  all  the  variables,  not  only  will  we  have  a richer  variety  of 
models  at  hand,  but  we  can  also  give  a better  interpretation  of  the 
table  under  consideration. 
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The  contents  of  this  thesis  will  be  outlined  in  this  chapter  to 
clarify  the  need  for  specialized  models  for  ordinal  data.  In  Chapter 
Two,  we  will  review  models  for  two-way  cross-classifications.  We  will 
discuss  models  where  both  variables  are  treated  as  nominal  in  scale, 
then  models  where  only  one  of  the  categorical  variables  in  the  table  is 
ordinal,  and  finally  models  for  tables  that  contain  two  ordered 
categorical  variables.  The  presentation  of  these  models  is  in  a similar 
vein  to  that  of  Goodman  (1979),  except  that  Goodman  defined  models 
through  properties  of  odds  ratios,  while  we  define  them  in  the  form  of 
log  expected  frequencies.  Our  approach  helps  to  clarify  how  the  new 
models  relate  to  the  standard  loglinear  models.  We  will  derive  the  like- 
lihood equations  for  these  models,  and  describe  iterative  methods  for 
finding  the  maximum  likelihood  (ML)  estimates  of  the  model  parameters. 
Available  procedures  for  finding  those  estimates  will  be  reviewed. 

Finally,  an  example  will  be  presented. 

Chapter  Three  consists  of  generalizations  of  the  models  introduced 
by  Goodman  (1979)  and  estimation  procedures  for  two-way  tables,  so  that 
they  can  be  applied  to  multi-dimensional  tables.  Primary  emphasis  will 
be  for  three-way  cross-classifications.  Interpretation  of  models  and 
examples  will  be  given  special  emphasis. 

Chapter  Four  concerns  testing  for  independence  in  two-way  tables. 

Most  often,  the  Pearson  chi-square  and  the  likelihood-ratio  chi-square 
are  used  to  test  for  independence  of  the  row  and  column  variables  in  a 
two-way  cross-classification.  However,  when  the  row  and  column  variables 
both  have  ordered  categories,  Haberman  (1981)  employed  both  the  canonical 
correlation  analysis  of  Fisher  (1940)  and  a model  described  in  Chapter  Two. 
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Instead,  we  propose  test  statistics  based  on  a simpler  model  than  the 
one  considered  by  Haberman.  We  will  review  briefly  the  conventional 
chi-square  statistics  (Pearson  and  likelihood-ratio)  and  the  statistics 
proposed  by  Haberman.  Then  we  describe  our  proposed  statistics  with 
their  properties.  Finally,  we  conduct  a simulation  study  to  investigate 
the  performance  of  our  proposed  statistics,  comparing  them  with  the  more 
familiar  Pearson  and  likelihood-ratio  chi-square  tests  for  independence. 

In  the  logl inear  models  for  tables  that  contain  ordinal  categorical 
variables  we  assign  scores  to  the  ordered  categories.  A limitation  of 
these  models  is  that  different  researchers  might  assign  different  scores 
to  the  ordered  categories.  To  overcome  this  drawback,  we  will  also 
discuss  models  in  which  the  scores  are  parameters  to  be  estimated. 
However,  the  estimated  scores  for  these  models  need  not  be  monotonic. 

It  is  desirable  to  have  the  estimated  scores  monotonic  if  one  can  assume 
the  true  scores  to  be  monotonic,  especially  when  an  apparent  natural 
order  applies  to  all,  or  some,  of  the  categorical  variables  in  a con- 
tingency table.  In  the  final  chapter  of  this  thesis,  we  will  propose 
procedures,  for  two-way  cross-classifications,  which  produce  monotonicity 
of  the  estimated  parameter  scores. 


CHAPTER  TWO 

ASSOCIATION  MODELS  FOR 
TWO-WAY  TABLES 

In  this  chapter  we  review  models  proposed  by  Goodman  (1979)  for 

cross-classifications  which  contain  ordered  categorical  variables. 

Although  the  contents  of  this  chapter  appear  in  the  literature,  they 

are  reviewed  for  the  convenience  of  the  reader.  In  section  1,  models 

defined  in  the  form  of  log  expected  frequencies  will  be  presented  and 

interpreted.  Maximum  likelihood  estimation  of  the  model  parameters  is 

discussed  in  section  2.  Section  3 describes  Newton-Raphson  procedures 

for  solving  the  likelihood  equations.  Also  a brief  review  of  available 

procedures  for  finding  the  ML  estimates  of  the  model  parameters  will  be 

given.  Finally,  an  example  will  be  used  for  illustration  in  section  4. 

Suppose  we  have  a contingency  table  with  r rows  and  c columns. 

The  number  of  observations  in  cell  (i,  j),  denoted  by  n..  (l^i<r,  l<j<c), 

^ J 

is  a realization  of  a stochastic  variable  with  expectation  m...  Hence- 

■ J 

forth  the  symbol  df  will  denote  the  residual  degrees  of  freedom.  Now, 
df  is  equal  to  the  difference  between  the  number  of  cells  in  the  table 
and  the  number  of  identifiable  parameters  under  the  fitted  model.  As  a 
goodness  of  fit  statistic  for  the  stated  model,  we  will  use  the  likelihood- 
ratio  statistic  G2,  defined  as 

G2  = 2ln..  log(nij./m.j ), 
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where  denotes  the  estimated  expected  frequency  in  cell  (i,j)  under 

the  stated  model.  A linear  model  in  the  natural  logarithms  of  the  cell 

frequencies  is  referred  to  as  a log! inear  model . Under  the  assumption  that 

o 

a certain  model  holds  the  statistic  G is  asymptotically  chi-square 

distributed  with  df  degrees  of  freedom. 

2 

The  use  of  G also  enables  us  to  test  the  validity  of  a certain 
logl  inear  model  M2  under  the  assumption  that  the  logl  inear  model  M-j 
holds.  The  model  M-j  should  be  a submodel  of  M2-  Denote  the  values  of 
G2  for  M-j  and  M2  by  G2(M1)  and  G2(M2),  respectively.  Then  under  the 
assumption  that  model  M-j  holds 

g2(m1|m2)  = g2(m1)  - g2(m2) 

is  asymptotical ly  chi-square  distributed  with  degrees  of  freedom 

d^l  ■ df2>  where  df^  and  df2  are  the  number  of  degrees  of  freedom  for 

the  residual  under  models  M-j  and  M2  respectively  (Bishop  et  al.,  1975,  p.525). 
2 

G (M1 |M2)  is  referred  to  as  the  conditional  test  statistic,  and  can  be 
interpreted  as  a measure  of  the  distance  between  models  M-j  and  m2> 

For  the  rxc  table  we  rarely  expect  the  independence  model 

log  m..  = y + A*  + \\  (ZA*  = = 0)  (2.1) 


for  which  df  = (r  - 1 ) (c  - 1),  to  give  an  adequate  fit.  If  model  (2.1) 

holds,  then  the  row  variable  (X)  and  column  variable  (Y)  are  said  to  be 

independent.  However,  in  the  scheme  of  the  nominal  hierarchical  models, 

the  model  of  next  greater  complexity  is  log  m^  = y + a*  + A^  + A^. 

The  model  has  (r  - l)(c  - 1)  additional  independent  parameters 
XY 

{^i j)  describing  the  association  between  the  categorical  variables 
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X and  Y.  The  model  has  as  many  independent  parameters  as  the  table  has 
cells  (i.e.  df=0)andis  referred  to  as  the  "saturated"  model.  Thus, 
among  the  usual  hierarchical  logl inear  models,  no  nontrivial  models 
exist.  Both  the  independence  and  saturated  model  treat  the  categorical 
variables  as  if  they  were  in  nominal  scale.  Thus,  the  above  models  are 
invariant  with  respect  to  changes  in  order  in  any  set  of  categories,  in 
the  sense  that  the  parameter  estimates  and  chi-square  statistics  remain 
unchanged  if  the  rows  or  columns  are  permuted.  If  one  or  both  variables 
are  ordinal,  however,  unsaturated  models  exist  that  are  more  complex  and 
usually  more  realistic  than  the  independence  model.  We  shall  discuss 
these  below,  first  for  the  case  of  a two-way  table  with  one  ordinal 
variable,  then  for  the  case  of  a two-way  table  with  both  categorical 
variables  ordered. 

2.1 . Association  Models 

Suppose  that  the  column  variable  (say  Y)  of  a two-dimensional  table 

is  ordinal  and  that  we  choose  a priori  a set  of  scores  {v.},  where 

J 

v1<v2<..<vc.  A model  that  is  more  complex  than  the  independence  model 
can  be  expressed  as  follows: 

log  miJ  = y + + (vj--v)Bi.  (2.2) 

Without  loss  of  generality  we  may  take  = zx\  = zb.  = 0.  Otherwise, 

J 

if  EB.,-  = B for  instance  we  could  redefine  "b.  = $.  - B/r,  ~X*.  = x\  + (v.-"v)B/r 

J J J 

in  order  to  satisfy  the  given  constraints.  The  model  has  df  = (r  - 1) 

(c  - 2),  and  is  unsaturated  when  r>l  and  c>2. 

Note  the  association  term  (2-factor  interaction)  (v.-V)B-  depends 

J 

on  the  ordinal  variable  Y through  the  assigned  column  scores  { v . } . 

J 
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The  independence  model  is  a special  case  of  model  (2.2)  with  all  6i  = 0. 

Within  a particular  row  (fixed  i),  the  deviation  of  log  m. . from  the 

■ J 

corresponding  quantity  in  the  independence  model  (2.1)  is  a linear  function 
of  the  ordinal  variable  (score)  with  slope  $...  If  $..>0  (6^0),  we 
expect  more  observations  above  (below)  the  mean  of  Y in  row  i than  if 
the  variables  were  independent.  For  any  pair  of  rows  i and  i ' , 
log(mij./mi , ^ ) = (X*  - X*. ) + (vj-7)(e.  - B.,). 

That  is,  on  the  log  scale  the  ratios  of  proportions  for  any  two  rows  is 

a linear  function  of  the  ordinal  variable  with  slope  c2(i,  i‘)  = 

i i 1 

(B.j  - B.j  i ) • Let  0.. , denote  the  odds-ratio  of  rows  i,  i'  and  columns 
j . J ' • 


• - , m. ./m. , • 
e11,  = -1J  i J 


m.  .m. , . , 
= AJ  i'J1 
m. . ,m. . . 
1J  i ' J 


ii  * 

Then  log  0j ^ ■ = (B^  - B1-i)(vj  - v^.,).  Hence  the  log  odds  ratio  is 

proportional  to  the  difference  between  column  scores.  For  integer  scores 

(Vj  = j}  the  odds  ratio  is  constant  in  rows  i,  i'  and  equals  exp  (B^  - B^- . ) 

for  all  pairs  of  adjacent  columns. 

Unlike  standard  loglinear  models,  model  (2.2)  is  invariant  only 

under  changes  of  rows.  For  integer  scores  {v.  = j}  Goodman  refers  to 

J 

the  model  as  the  "row-effect"  model,  and  we  refer  to  it  as  the  R model . 

Suppose  now  that  the  row  variable  is  ordinal  and  the  column  variable 

is  nominal,  and  suppose  that  we  choose  a priori  a set  of  scores  {u^} 
where  u-|<u2<. . .<ur«  A parallel  model  to  the  R model  is 

log  m..  = p + X*  + X^  + (u.-ujBj. 

Here  EX^  = = ZBj  = 0 and  df  = (r  - 2)(c  -1). 


(2.3) 
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All  results  obtained  for  the  R model  can  be  translated  directly  for 
model  (2.3).  For  integer  scores  (ui  = i } , Goodman  refers  to  it  as  the 
"column  effect"  model  and  we  refer  to  it  as  the  C model. 

Now  let  both  the  row  variable  X and  the  column  variable  Y of  the 
table  be  ordinal.  We  suppose  that  we  choose  a priori  two  sets  of 
scores  {u^}  and  { v^ } to  assign  respectively  to  the  rows  and  columns, 
where  u-j<U2<...<ur  and  Vi<v2<. . .<vc.  A parsimonious  logl  inear  model, 
which  utilizes  the  ordinal  information  and  has  only  one  more  parameter 
than  the  independence  model , is 

log  mij.  = y + X*  + X^  + (u^u) (Vj-v)B.  (2.4) 

The  model  has  residual  df  = rc  - r - c and  is  unsaturated  for  all  but 
2x2  tables.  The  independence  model  is  a special  case  with  8=0. 

If  6>0,  there  are  more  large  X-large  Y and  small  X-small  Y observations 
than  would  be  expected  if  X and  Y were  independent.  The  deviation  of 
log  m^  . from  the  corresponding  quantity  in  the  independence  model  is  a 
linear  function  of  the  row  scores,  for  fixed  level  j of  Y,  with  slope 
(Vj-v)B,  and  a linear  function  of  the  column  scores,  for  fixed  level 
i of  X,  with  slope  (u.j-u)B.  Thus,  we  refer  to  the  model  as  the  li near- 
by-linear  association  model.  For  any  pair  or  rows  i and  i‘  and  any  pair 
of  columns  j and  j' 

leg  e]’!  = e(u1-ui.)(vj-vj.). 

That  is,  the  farther  the  rows  and  columns  are  apart,  the  greater  in 
absolute  value  the  log  odds  ratio  becomes.  For  the  scores  {u^  = i}  and 
(Vj  = j}  the  odds  ratio  is  constant  and  is  equal  to  exp  (s)  for  the 
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(r  - 1 ) (c  - 1 ) 2 x 2 subtables  of  the  r x c table  formed  from  adjacent 
rows  and  adjacent  columns.  Due  to  the  constancy  of  the  odds  ratio 
Goodman  refers  to  it  as  the  "uniform  association"  model. 

Model  (2.4)  is  not  invariant  to  permutations  of  rows  and  columns. 

The  model  is  appropriate  if  the  li near-by-linear  association  property 
holds  for  the  set  of  scores  chosen  by  the  user.  Otherwise,  one  can 
alternatively  treat  the  scores  {u . } and  {v^}  as  parameters  {y i>  and 
{vj}  to  be  estimated  for  which  1 inear-by-1 inear  association  is  best 
approximated.  Such  a model  can  be  expressed  as  follows: 

log  m..  = y + X*  + + y.jVj3*  (2.5) 

Here  we  may  take  ZX.j  = ZXj  = Zy.j  = ZVj  = 0.  Further  arbitrariness  in 

the  score  parameters  is  removed  by  scaling  them  and  setting,  without 

loss  of  generality,  Zy2  = Zv2  = 1.  Otherwise  if  Zy2  = y2  and  Zv2  = v2 

J i j 

for  instance,  we  could  redefine  6 =6yv,  ^ = y./y,  v-  = v./v  in  order 

to  satisfy  the  additional  constraint  above. 

Unlike  previous  models,  model  (2.5)  is  multiplicative  rather  than 

linear  in  the  parameters.  It  also  differs  from  the  1 inear-by-1 inear 

association  model  in  the  sense  that  it  is  invariant  under  changes  of 

rows  and  columns.  However,  if  the  estimated  scores  are  monotonic,  one 

can  interpret  the  model  through  odds  ratios  merely  by  replacing  the 

fixed  scores  {u.}  and  { v . } by  the  parameter  scores  {y.}  and  {v.} 

J 1 J 

respectively.  The  relative  distances  between  the  rows  and  columns  in 
order  for  the  odds  ratio  to  be  constant  per  unit  distance  are  indicated 
by  the  distances  between  the  estimated  scores.  If  only  one  set  of 
scores  is  in  order,  say  the  column  scores  {v.},  then  the  row  scores  {y.} 
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can  be  interpreted  as  row  effects  as  in  the  R model  (2.3).  That  is,  in 
that  case  the  model  can  be  viewed  as  the  R model  by  setting  By..  = 3- 
in  equation  (2.5).  Now,  if  both  sets  of  scores  are  nonmonotonic,  then  the 
row  scores  and  column  scores  can  be  interpreted  as  row  and  column 
effects  respectively.  In  that  case,  the  association  between  the 
variables  is  positive  in  some  locations  and  negative  in  others.  In 
Chapter  Five  we  will  propose  a procedure  for  estimating  the  score 
parameters  subject  to  the  constraint  that  they  are  monotonic.  This 
model  has  been  discussed  by  Andersen  (1980)  and  Goodman  (1979),  who 
refers  to  it  as  "Model  II,"  and  we  refer  to  it  as  the  RC  model . One 
final  note  about  the  RC  model  is  that  the  conditional  likelihood-ratio 
goodness  of  fit  statistic  G^(I|RC)  doesn't  have  an  approximate  chi-square 
distribution  when  r>2  and  c>2  (Haberman,  1981).  The  statistic  G2 ( I | RC ) 
tests  for  independence  given  that  the  RC  model  holds. 

The  final  model  discussed  is  a generalization  of  all  models  dis- 
cussed above,  except  for  the  RC  model.  It  can  be  expressed  as  follows: 

log  mij  = y + X*  + x]  + (u.-u) (v-v)bXY  + (v^-vjB*  + (u.-u)Bj,  (2.6) 

where  ZXX  = IXY  = 0 and  B*  = 3*  = b]  = BY  = 0. 

The  model  has  residual  df  = (r  - 2)(c  - 2),  and  is  unsaturated  whenever 
r>2  and  c>2.  To  clarify  the  constraints  on  the  B-parameters,  we  need  to 
look  at  the  local  log  odds  ratios.  Suppose  that  we  use  integer  scores 
for  both  sets  of  scores.  For  all  pairs  of  adjacent  rows  and  adjacent 
columns  we  have 


^ - sXY 


<8i+X>  + (bJ+1-»J)  • eXY  ♦ B*  ♦ Bj. 
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Thus  the  local  log  odds  ratio  is  a linear  function  of  differences  be- 
tween the  relevant  row  effects,  $?,  and  differences  between  column 

C 

effects,  Bj-  The  usual  kinds  of  conditions  included  in  additive  models 
r-1  R c-1  c 

are  z 3-  = £ 3,-  = 0,  and  this  constraint  is  equivalent  to  3?  = rX 

i=l  1 j=l  J 1 r 

and  3^  = 3^.  Without  loss  of  generality  we  may  take  3X  = 0 and  3^  = 0. 

If  3*  * 0 for  instance,  we  could  redefine  ]3X  = 3X  - 3*  and 
— Y Y X — 

= + 6l^vj“v^  in  order  t0  satisfy  the  constraints. 

Goodman  refers  to  model  (2.6)  as  "Model  I,"  and  we  refer  to  it  as 
the  R+C  model . For  further  details  on  interpretations  on  these  and 
similar  models  see  Haberman  (1974a),  Simon  (1974),  Goodman  (1979,  1981) 
and  Agresti  (1983). 


Summary  of  Models 

Page 

Independence  model 

(2.1) 

5 

R model 

(2.2) 

6 

C model 

(2.3) 

7 

Linear-by-1 inear  association  model 

(2.4) 

8 

RC  model 

(2.5) 

9 

R+C  model 

(2.6) 

10 

2.2.  Maximum  Likelihood  Parameter  Estimation 

The  parameters  in  log! inear  models  can  be  estimated  by  the  method 
of  maximum  likelihood.  To  write  down  the  likelihood  equations, 
assumptions  about  the  sampling  distribution  must  be  made.  Three  possible 
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schemes  will  be  discussed:  independent  Poisson,  multinomial  and  product 

multinomial.  We  will  present  these  for  two-dimensional  contingency  tables. 
If  the  observations  n..(l<i^r,  l<j<c)  are  realizations  of 

' J 

independent  Poisson  random  variables  with  parameters  m..,  the  likelihood 

i J 

can  be  written  as 


it  e 
i,j 


-m. . m. . 
U _U 


n. .! 


In  the  case  of  the  multinomial  sampling,  the  sample  size  N is 

fixed  and  the  vector  (n^ nrc)  is  a multinomial  random  variable, 

' ' 

with  probability  vector  Vl_  with  N = Em.  . = Zn...  The 

N ’**•’  N 


with  N = Em.  . = Zn . 

1J  ij 


likelihood  L..  becomes 
M 


Lm  = N!  it 


m. 


21 


n . . 
i,J  iJ 


n . . 
U 


Under  product  multinomial  sampling  either  the  row  totals  or  column 

totals  are  fixed.  With  fixed  row  totals  n.+  = zm.  . = zn..,  it  is 

1 j 1J  j 1J 

assumed  that  the  vectors  (n^,...,n^c),...,(nr^,...,n  ) are  independently 

multinomial  distributed  with  parameter  vectors 


The  likelihood  Lp^  becomes 


mil 


V 


m. 

ic 

'v 


for  1 <i <r . 


lPM  ^ni  + ! 


m.  . 
_LI 


V 


n . 

ij 


ij 


Now,  given  a dataset  {n^.}  and  model,  maximum  likelihood  estimates  under 
each  of  the  three  sampling  schemes  can  be  obtained.  Haberman  (1974b,  p.  21) 
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shows  that  the  ML  estimates  exist  if  all  n..>0  l<i<r,  1 <j . Denote  the 

■ J 

ML  estimates  of  the  independent  Poisson,  multinomial,  and  product 
multinomial  by  m^,  m^.,  and  m^.  respectively.  By  rewriting  LM  and  LpM 
as 


LM 


N!e 


-m. . n . . 
e 


ij  niJ! 


N 'e^ 


lpm  " 71 


n . ip 
i+'e 


v 


i w 


( 

-m . . n . . 

ni+] 

/1+e 

. n.  .i 

J U' 

i ni+ni+ 

t 

V 


it  follows  (Bishop  et  al . , Thorem  13.4.1)  that 


AP 

if 

„ . . 

mij 

= m.  . 
ij 

In,  = N 
lj  IJ 

(2.7) 

. ~P 
and  m. . 
ij 

- £pm 
- m. . 

ij 

if 

Zitl . = n . 
j 1J  1 + 

Assume  one  of  the  sampling  schemes,  say  Poisson  sampling.  Then  the 

1 ogl i kel i hood  function  K (m)  is 

P “ 

Kp(m)  = Log  Lp  = En.j  log  m^.-Zm^- z log  n^.!.  (2.8) 


To  derive  the  likelihood  equations,  we  consider  the  R+C  model  first, 
since  it  is  the  most  general  of  the  loglinear  models  already  discussed. 
Then  we  consider  the  RC  model.  Substituting  (2.6)  into  (2.8)  we  obtain 
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Kp(— ) " + ^i+^i  + ^+.\ Y + eXYE(ui-u)(vj.-v)nij  + z(vj-V)nij6^ 

+ 2(vj'v)ni jBj  - Eexp(y+AiX+Aj+BXY(ui-u)(Vj-v)  + (vj-v)BX 
+ (ui ~u )Bj ) - Z log  n 

The  likelihood  equations  are  obtained  from  the  following  partial 
derivatives: 


3Kp(m) 


ni+  - Zexp(y+AX+AY+BXY(uru)(v.-v)  + bJ(v.-v)  + BY(u,-u)) 
i j 1 j * j j i 


X,  -x  . JT,  -> 


ni+  " mi+ 


1 <i<r. 


(2.9a) 


By  symmetry  we  have 


9K  (m) 

— Y — = n . -•  * in, 


3A 


+J  +j! 


J 


(2.9b) 


9Kp(m) 


5(¥j-v,nu 


2(vrv)  exp  [log  mj  .]  = E(v.-v)(n4  .-m.  .) 


U U 


l£i£r. 


(2.9c) 


By  symmetry  we  have 


3K  (m) 

— Y = 2(uj"u)(nii_mii); 
3B-.  i 1J 


(2.9d) 


3K  (m) 


3B‘ 


XY  “ 2(ui-u)(vj.-v)nij.  - Z(ui-u)(vj.-v)  exp  [log  m^] 


.jI(ui-u)(vrv)(nirmij) 


l.J 


(2.9e) 
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We  obtain  the  likelihood  equations  by  setting  (2.9,  a-e)  equal  to  zero. 
Note  that  equation  (2.9e)  is  implied  by  either  (2.9c)  or  (2.9d). 
Therefore  the  ML  estimates  for  the  parameters  of  the  R+C  model  are 
given  as  the  solutions  to 


A 


mi+  ~ ni+ 

1 <i <r ; 

(2.10a) 

m+j  = n+j 

l^j<c ; 

(2.10b) 

= Z ( v .-v)n . . 
j J ’j 

1 <i^r ; 

(2.10c) 

Z(u,-u)m  = 
i J 

l^j^c. 

(2.1 Od) 

Subsets  of  the  above  constraints  are  satisfied  by  the  R model  and 
C model.  For  instance,  the  ML  estimates  {m..,  Isi^r,  l<j<c}  of  the 
expected  frequencies  under  the  R model  satisfy  the  set  of  equations  (2.10a) 
through  (2.10c).  However,  for  the  1 inear-by-1 inear  association 
model,  the  ML  estimates  (m..,  l<i<r,  l<j<c},  in  addition  to  equations 
(2.10a)  and  (2.10b),  should  satisfy  the  following  equation  obtained  by 
setting  equation  (2.9e)  equal  to  zero: 


Z(u.-u)(v.-v)m ..  = Z(u.-u)(v,-v)n. .. 
i » j J J i*j  1 J 1J 


(2.10e) 


It  is  easily  seen  that  the  ML  estimates  {in..;  l<i<r,  1 < jsc > of  the 
expected  frequencies  under  the  RC  model,  in  addition  to  equations  (2.10a) 
and  (2.10b),  should  satisfy  the  following  equations  obtained  by 
differentiating  the  loglikel ihood  function  and  setting  the  partial 
derivatives  equal  to  zero: 


Zv -m. . = Zv.n . . 

jJiJ 

Zu.m. . = Zp.n . . . 
i i U j l 


(2.1  Of ) 
(2.10g) 
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The  Independence  model  (2.1)  satisfies  (2.10a)  and  (2.10b). 

That  is,  the  observed  and  fitted  frequencies  have  the  same  row 
marginals  and  the  same  column  marginals.  Since  model  (2.1)  is  a 
special  case  of  all  models  discussed,  the  expected  frequencies  under 
all  models  share  the  same  above  property.  For  instance,  in  addition 
to  that,  the  observed  and  fitted  row  means  under  the  R model  (2.2) 
are  equal.  From  equation  (2.10e)  we  see  that  the  observed  and  fitted 
moments  of  the  form  E(XY)  are  equal.  Hence,  for  the  1 inear-by-1 inear 
association  model,  we  see  from  equations  (2.10a,b,e)  that  the  observed 
correlation  equals  the  fitted  correlation.  For  higher  dimensional 
tables  (next  chapter)  we  will  refer  to  moments  of  the  form  E(XYZ)  as 
moments  of  third  order  or  third  moments. 

We  remark  from  (2.7)  and  (2.10)  that  the  maximum  likelihood 
estimates  of  the  expected  frequencies  under  all  the  models  discussed 
previously  are  the  same  under  each  of  the  three  sampling  schemes 
(Poisson,  multinomial,  and  product  multinomial). 

2.3.  Newton-Raphson  Procedures 

The  solution  of  the  equations  (2.10)  is  a numerical  analysis 
problem.  In  this  section  we  discuss  iterative  methods  designed  for 
solving  the  set  of  equations  (2.10). 

Suppose  we  want  to  solve  a set  of  p nonlinear  equations  with  p 
unknown  variables 

i.e.:  f . (x ) = 0 l<i<p  where  x = (x1 ,x2> . . . ,xp) 

or  f(x)  = ( f 1 (xj , f2(x) f (x))  = 0. 
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Let  J = { J j , l<i,  j<p)  denote  the  Jacobian  matrix  of  f (xj  and  H = 

{^ij»  »j^p}  its  inverse,  assuming  J to  be  nonsingular  where 

,(x)_  3f1  (*> 
ij  3x. 

J 

The  classical  Newton  method  is  obtained  by  linearizing  f.  If  we  assume 
x = a is  a zero  for  f,  that  x^  is  an  approximation  to  a,  and  that  f is 
differentiable  at  x = x^,  then  to  a first  order  approximation 

0 = f(a)sf(xQ)  + J^Ma  - Xq). 

If  J (jCq)  is  nonsingular,  then  the  equation  ffxg)  + J^)^  - Xq)  = 0 
can  be  solved  for  x^ : 

2L]  ~ 2Lq  ” ^ f ( Xq  ) 

and  x ^ can  be  taken  as  a closer  approximation  to  a the  zero  of  f. 

Therefore,  the  Newton-Raphson  iterative  method  for  solving  the  above 
system  of  equations  is  given  by 

2Li+l  = 2L-j  - H(x.)f(x.)  i = 0,1,2,...  (2.11) 

where  x.  is  an  approximation  to  the  zero  a at  the  ith  iteration. 

We  note  that  for  the  case  of  a function  of  one  variable,  a special 
case  of  (2.11)  referred  to  as  the  unidimensional  Newton-Raphson  method,  can  be 
expressed  as 

xi+1  = X.  - f(x.)/f,(x.)  , (2.12) 

where  f1  is  the  derivative  of  f.  The  unidimensional  Newton-Raphson 
iterative  method  (2.12)  can  be  applied  to  solve  a set  of  p equations 
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with  p unknown  variables.  To  do  so,  let  ^ be  a chosen  starting  value, 

then  consider  ^ (x ) as  a function  of  x-j  only,  with  the  remaining  variables 

fixed  at  values  ^20*^30’ ** **p0*  ^ solve  f-|(x-|)  = 0 to  obtain  an 

approximate  zero  x^.  Then  consider  f 2 (ii)  as  a function  of  x^  only,  with  the 

remaining  variables  fixed  at  values  x^.x^.^.x  Q.  We  solve  f2(x2)  = 0 

to  obtain  an  approximate  zero  x21 . We  continue  this  process  until 

convergence  is  reached,  in  the  sense  that  ||f(x)||  = max  |f.(x)|<e, 

l<i<p  1 

for  tolerance  e chosen  by  the  user. 

Before  illustrating  the  application  of  the  iterative  methods  (2.11) 
and  (2.12),  we  should  mention  their  properties.  First,  method  (2.12) 
does  not  require  inversion  of  matrices  (i.e.;  H,  the  inverse  of  the 
Jacobian  J).  However,  both  methods  were  used  for  many  examples  and  it 
was  found  that  method  (2.12)  required  many  more  iterations.  Finally, 
in  our  own  case  the  function  f consists  of  the  partial  derivatives  of 
the  loglikelihood  function;  therefore,  J is  the  information  matrix  and 
-H  will  be  the  asymptotic  variance-covariance  matrix.  Hence  the 
advantage  of  method  (2.11)  is  that  it  provides  an  estimate  of  the 
variance-covariance  matrix  of  the  parameter-estimates,  through  -H  = 
r lim  - H(xf)). 

As  will  be  seen  later,  the  form  of  the  Jacobian  is  simple  when 
the  model  is  logl inear.  Hence,  for  ease  of  presentation  we  will 
illustrate  the  application  of  method  (2.11)  for  loglinear  models  and 
method  (2.12)  for  the  multiplicative  model. 

H - (n-j  -j , . . . ,nrc)  be  the  rc-vector  of  observed  counts,  where 
the  {n^}  are  realizations  of  independent  Poisson  random  variables  with 
parameters  (m^}.  Let  m = (m^  ,. . . ,mrc)  be  the  corresponding  rc-vector 
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of  the  expected  frequencies  under  the  fitted  model.  Let  Np  and  Mp 
denote  respectively  the  diagonal  matrices  with  observed  and  expected 
counts  on  the  main  diagonal.  Let  be  the  p-vector  of  parameters. 
Our  model  can  be  expressed  as 


Log  (m)  = Xb 

where  X is  the  design  matrix  which  has  form  depending  on  the  model. 

It  can  be  seen  from  equations  (2.9)  that  the  function  f for  logl inear 
models  is  (n-m)'X  and  the  Jacobian  J is  as  follows 


J ■ ft  - 


where  m = exp  [XbJ. 


Now 


xnmn xiPmri 

o 

• 

• 

CM 

,e" 

* ' 

xllx12  xlp 

Xrc, lmrc  xrc,pmrc 

o 

• 

• 

3 

-s 

o 

xic,r-*-xic,p. 

Therefore  J = -X'MpX.  Using  equation  (2.11)  we  obtain 

bi+1  = b.  + (X'M^Xf  Vfn-m.)  i = 0,1,2,...  (2.13) 

where  m^  = exp  [Xb^.]  and  Mp  is  the  diagonal  matrix  with  the  estimated 

expected  frequencies  m.  for  iteration  i on  the  diagonal.  Like  any 
iterative  method,  a starting  value  is  needed,  which  may  be  taken  as 
follows: 


bp  = (X*NdX)_1X,Nd  l.og(n) 


mp  = exp  [Xbp] . 


(2.14) 
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Note  if  some  n^.  0,  one  could  add  1/2  to  each  cell  count  and  replace 

H in  (2-14)  by  in  + \ 1_  where  j_  is  a rc-vector  containing  l's.  Now 

iterate  with  equation  (2.13)  until  maxi  m. . ^t+1  ^ 

l^r  1J  U 1 ’ 


The  tolerance  e is  up  to  the  user  (e.g.:  e = .001).  To  illustrate 
the  above  method  consider  the  R model 

log  mij.  = u + \\  + A^  + (Vj-vJg.  (EA^  = EA^  = Eg.  = 0). 

Here,  if  T denotes  the  transpose  of  a vector,  we  have 

x X Y Y T 

k=  (y  ••  »Ar_i  »A-|  »•••  »AC_-|  »3i  »6r_i ) > a (2r+c-l ) -vector , 

— " ^nll’****nic’n21’***  ,n2c’  * ’ ‘ ,nrl  ’ * * ’nrc^’  a rc-vector, 
m = (m^  , . . . ,m^c,m2-| ,. . . ,m2c ny.j , . . *mrc)T,  a rc-vector,  and 


X = 

(rcx2r+c-l ) 


1 1 . 


1 1 
1 0 


1 0 
1 -1 


1 -1 


0 1 
. 0 

. 0 

0 -1 


0 


0 

-1 


1 

0 


1 
0 
• 

0 

1 -1 


-1 


1 

0 


0 

-1  -1 


0 

0 

1 

■1 


1 

-1 

0 


( V-i  —v ) 


(V-V) 

c0 


0 

0 


0 

0 


0_ 

-(vrv) 


"(vc-v) 


0 

0 


(v 


-V) 


(v  -v) 

-(Vi-V) 


-<vc-v) 
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Consider  now  the  RC  model  (2.5).  We  use  this  model  to  demonstrate 
the  unidimensional  Newton-Raphson  procedure.  If  we  exponentiate 
equation  (2. 5)and  reparametrize  it,  we  obtain  the  form  used  by  Goodman 
(1979). 


If  we  treat  equations  (2.10a),  (2.10b),  (2. 1 Of ) , and  (2.10g)  as 
pertaining  to  the  parameters  a,,  b.,  p.  and  a.  respectively,  then  we 
can  use  the  unidimensional  Newton-Raphson  method  to  find  the  ML  estimates 
of  the  parameters.  This  can  be  done  by  rewriting  equations  (2.10)  as 
follows: 


(2.10a)«-*-f.| .. (a. ) - 0, 


(2.10bW2j(bj)  = 0, 


(2.10f)~f3i(p.)  = 0, 


To  solve  the  equations  above,  we  proceed  as  follows: 


fl  i (ai ) " ni+“mi+  " = ai  - ciiV  where  cli  = mi+/ai  * 

Hence  f-j (a^ ) = -c^,  and  at  iteration  (t+lj.aj^  can  be  replaced  by 


(t+1)  , 

a.  , where 
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.(t+D  _ (ti  fii(ai  1 _ (t)  . ni+'mi+  ft)  , 

1 1 1 1 i+  i+ 


(2.13a) 


By  symmetry  b!^  can  be  replaced  by  b(t+1),  where 
J J 

,(t+l)  k(t)  , 

bj  = bj  W 


(2.13b) 


Now  for  equation  (2.1  Of)  we  have 


f31  (p1  > = £oj(nU-m1j)  - Zoj(nij-aib/'°J)-  a"d  f3i<Pi>  “ - 


EaZm.  .. 
J U- 


Thus  at  iteration  (t+1),  can  be  replaced  by  y|t+^,  where 


(t), 


,(t+i)  . (t)  f zqj  i(ni.i-mi.i) 


v 2TtT 

Zaj  mij 


At) 


By  symmetry  for  a.,  we  have 

J 


(2.13c) 


v 


v (t+1 ) / » 

(t+D  ...(0  . Sbi 

‘ vj  + r.2(t+l), 


£p 


mi.j 


(2.13d) 


Now  the  algorithm  given  by  Goodman  (1979)  goes  as  follows: 

Take  as  starting  values  aj0)  = bj0)  = 1,  yj0)  = i and  vj0)  = j so  that 

P-j  ^ = i - (r+l)/2,  and  oj^  = j - (c+l)/2.  At  stage  t of  the 

iterative  process,  let  the  estimated  parameters  be  denoted  by  a!^, 

bj  yi  anc*  At  stage  (t+1),  a given  parameter,  say  a j , 

can  be  replaced  by  ajt+1 ^ according  to  equation  (2.13a).  The  estimated 

expected  counts  at  that  stage  (t+D  are  m.  . = a!t+1^b!t^  exp  lo^a^h 

ij  i j i j 

Then  bj  \ pj  ^ and  aj  ^ can  be  replaced  in  turn  by  their  next  values. 
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and  the  estimated  expected  counts  are  computed  each  time  one  of  the 
sets  is  replaced.  Repeat  the  above  until  satisfactory  convergence 
occurs. 

Clogg  (1982)  gives  a Fortran  program  for  finding  the  ML  estimates 
of  the  parameters  of  al 1 the  models  discussed,  using  the  unidimensional 
Newton-Raphson  method  (2.12).  We  have  written  a program  that  uses  SAS 
(PROC  MATRIX)  for  finding  the  ML  estimates  of  the  parameters  in  log- 
linear  models,  such  as  the  R model  (see  Appendix  2). 

A third  method  for  finding  the  ML  estimates,  referred  to  as  the 
iterative  scaling  method  for  loglinear  models,  follows  from  Theorem  1 
of  Darroch  and  Ratcliff  (1972).  At  each  stage  of  the  iterative  process, 
the  method  approximates  the  expected  counts  under  the  stated  model. 

Let  m!j)  be  the  approximation  for  m^.  and  the  tth  stage.  Then  a single 
cycle  has  three  steps.  If  we  apply  the  procedure  to  the  R model,  we 
have  for  t = 0,1,2,... 


(t+1 ) 
m. . ' 

1J 


= ni+/m 


(t) 

i+ 


m 


(t+2) 

ij 


= n+j/m 


(t+1) 

+j 


(t+3)  (t+2) 

mij  = mij 


svi^v;e2) 


1 -V  . 

J 


★ 

where  the  {vk>  are  a linear  rescaling  of  the  column  scores  {v.  > that 
* ^ 
satisfy  0<vk<l.  By  symmetry  the  same  method  can  be  applied  to  the  C 

model.  For  the  1 inear-by-1 inear  association  model,  all  the  {v*}  in  the 

K 

above  formulae  should  be  replaced  by  u^v^  and  the  summation  is  over  l 
and  k.  To  illustrate  the  application  of  previous  models,  we  give  an 
example  next. 
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2.4.  Example 

The  6x4  cross-classification  in  Table  2.1  gives  the  relationship 

between  the  ordinal  variable,  mental  health,  and  ordinal  variable, 

parent's  socio-economic  status,  studied  by  Goodman  (1979).  The  values 

in  parentheses  consist  of  ML  estimates  of  expected  counts  for  the 

uniform  association  model  (2.4),  assuming  one  of  the  standard  sampling 

schemes.  Goodman  fits  all  models  discussed  previously  using  integer 

scores.  The  results  are  displayed  in  Table  2.2.  We  can  see  from  the 

table,  except  for  the  independence  model,  the  association  models 

fit  the  data  very  well.  Since  the  uniform  association  model  (2.4)  is 

the  simplest  model  among  those  fitting  the  data,  it  serves  well  as  a 

description  of  the  table.  As  noted  earlier,  whenever  a model  is  a 

submodel  of  another  model  (i.e.  the  latter  is  a generalization  of  the 

former),  we  can  obtain  a conditional  breakdown  of  the  likelihood-ratio 
2 

statistic  G to  test  the  goodness  of  fit  of  a model,  given  that  its 
generalization  holds.  Results  of  this  type  are  given  in  Table  2.3. 

Since  the  uniform  association  model  is  a submodel  of  the  R+C  model,  we 
can  test  for  row  and  column  effects.  From  the  table  we  see  that 
G (R+C)  - G (U)  = 6.85  based  on  df  = 6 is  not  statistically  significant. 
This  confirms  that  the  uniform  association  model  suffices  to  describe 
the  data.  Assuming  that  the  uniform  association  model  holds,  we  may 
test  for  independence  using  the  conditional  likelihood  chi-square 
statistic  G2(I|U)  = G2(I)  - G2(U)  = 37.53,  based  on  df  = 1.  The  test 
strongly  suggests  that  the  two  variables  are  dependent. 

The  uniform  association  model  gives  a good  fit,  with  G2  = 9.89 
based  on  df  = 14.  We  note  that  by  adding  a single  parameter  to  the 
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independence  model,  we  gain  a large  reduction  in  the  G2  statistic. 

The  estimated  constant  value  of  the  odds  ratio  is  found  to  equal  1.095. 
That  means,  for  instance,  that  the  odds  of  having  moderate  symptom 
formation  instead  of  impaired  mental  health  is  1.095  times  larger  for 
people  having  parent's  SES  in  category  C rather  than  category  D. 

The  same  can  be  said  for  any  pair  of  adjacent  rows  and  adjacent  columns. 

Let  us  consider  certain  invariance  properties.  If  we  interchange, 
say,  column  A and  column  B of  the  6x4  table,  the  chi-squared  values  will 
remain  unchanged  for  models  (2.1),  (2.3),  and  (2.5).  Only  for  the 
independence  model  and  the  RC  model  will  the  chi-squared  values  be 
invariant  under  permutations  of  both  rows  and  columns. 

There  are  many  existing  computing  methods  that  produce  the  ML 
estimates  for  the  logl inear  models  discussed  previously.  We  shall 
now  discuss  some  of  them,  (i)  The  computer  package  GLIM  (Generalized 
Linear  Interactive  Modelling:  Baker  and  Nelder  (1978))  is  designed 

for  fitting  generalized  linear  models.  A generalized  linear  model  is 
a model  whose  dependent  variable  has  a distribution,  with  parameter  6, 
which  belongs  to  the  exponential  family,  a set  of  independent  variables 
xi»x2»...,xm  predicted  by  E(Y)  = Ea^.x^ , and  a linking  function  9 = f(EY) 
connecting  the  parameter  e of  the  distribution  with  the  {Y's}  of  the 
linear  model  (Nelder  and  Wedderburn  (1972)).  Since  the  distribution 
(Poisson)  considered  hereisa member  of  the  exponential  family, 
where  the  natural  parameter  (log  m)  is  a linear  function  of  the 
independent  variables,  the  logl inear  models  are  special  cases  of 
generalized  linear  models.  The  package  mentioned  is  applicable  to 
logl inear  models  discussed  previously.  It  uses  the  multidimensional 
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Newton-Raphson  procedure,  (ii)  The  computer  package  SPSSX  contains  a 
LOGLINEAR  procedure  which  performs  model  fitting,  hypothesis  testing 
and  parameter  estimation  for  any  model  that  has  categorical  variables 
as  its  major  components.  It  produces  ML  estimates  by  means  of  the 
Newton-Raphson  method  (2.11).  Output  includes  observed  and  expected 
cell  counts,  residuals  and  adjusted  residuals,  and  both  the  Pearson  and 
likelihood  ratio  chi-square  statistics.  One  can  request  printing  of  the 
design  matrix,  parameter  estimates,  standard  errors,  standardized 
parameter  estimates,  confidence  intervals  and  correlation  matrix  of  the 
estimated  parameters,  (iii)  Haberman  (1979)  gives  a Fortran  program 
(FREQ)  that  also  fits  logl inear  models  using  the  multidimensional 
Newton-Raphson  procedure,  (iv)  Duncan  and  McRae  (1979)  show  how  to 
fit  ordinal  logl inear  models  through  iterative  use  of  the  computer 
package  ECTA  of  Fay  and  Goodman  (1975),  which  is  based  on  iterative 
proportional  fitting. 

We  remark  that  all  the  available  computing  procedures  are  designed 
for  logl inear  models.  However,  they  can  be  applied  to  non-logl inear 
models.  For  instance,  they  can  be  applied  to  the  RC  model  in  the 
following  way: 

We  may  rewrite  the  model  as 

X y — — 

log  rru j = y + + y.\K,  where  y..  = y^. 

If  we  treat  the{v,}  as  known  scores  (e.g.  v.  = j - (c+l)/2),  then  we 
obtain  the  R model  which  can  be  fitted  via  one  of  the  above  procedures. 

A 

The  fitted  R model  produces  estimates  {y.}  of  the  y.  scores.  Then  we 

A 

treat  the  {y^}  as  fixed  to  obtain  the  C model  which  can  be  fitted  also 
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through  one  of  the  procedures  above.  The  fitted  C model  produces 
estimates  {v.}  of  the  v.  scores.  These  scores  will  be  used  to  fit  the 

J J 

R model  again.  The  procedure  fixes  one  of  the  set  of  scores  and 
estimates  the  other  set  until  adequate  convergence  is  reached. 
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Table  2.1.  Cross-classification  used  to  illustrate  the  linear-by- 
1 inear  association  model. 


Parent1 s 
socio-economic 
status 

Mental  Health  Status 

Well 

Mild  symptom 
formation 

Moderate  symptom 
formation 

Impaired 

A 

64 

94 

58 

46 

(65.3) 

(104.4) 

(50.2) 

(42.1) 

B 

57 

94 

54 

40 

(54.2) 

(94.9) 

(49.9) 

(45.9) 

C 

57 

105 

65 

60 

(55.9) 

(107.2) 

(61.7) 

(62.2) 

D 

72 

141 

77 

94 

(65.3) 

(137.0) 

(86.4) 

(95.3) 

E 

36 

97 

54 

78 

(39.0) 

(86.9) 

(61.8) 

(74.7) 

F 

21 

71 

54 

71 

(27.3) 

(68.8) 

(52.0) 

(68.8) 

Note:  The  parenthesized  values  are  the  estimated  expected  frequencies 

under  the  uniform  association  model. 


Table  2.2.  Association  models 

applied  to  Table  2.1. 

Association  models 

Degrees  of 

Likel i hood-ratio 

freedom 

chi-square 

(2.1)  Independence  model 

15 

47.42 

(2.4)  Uniform  association  model 

14 

9.89 

(2.2)  R-model 

12 

6.28 

(2.3)  C-model 

10 

6.83 

(2.6)  R+C  model 

8 

3.04 

(2.5)  Multiplicative  model 

8 

3.57 

Note:  The  parenthesized  numbers 

cussed  previously. 

are  the  corresponding 

models  dis- 
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Table  2.3.  Analysis  of 

association 

for  Table  2.1 . 

Effects  of  association 

Models  used 

Degrees  of  freedom 

Likelihood  ratio 
chi-square 

General  effect 

(2. 1 ) - (2. 4 ) 

15  - 14  = 1 

37.53 

Row  and  column  effects 

(2.4)-(2.6) 

14-8  =6 

6.85 

Other  effects 

(2.6) 

8 

3.04 

Total  effects 

(2.1) 

15 

47.42 

CHAPTER  THREE 

ASSOCIATION  MODELS  FOR  MULTIDIMENSIONAL 
CROSS-CLASSIFICATIONS  OF  ORDINAL  VARIABLES 

In  this  chapter  we  shall  extend  Goodman's  models  for  two-way 
tables,  so  as  to  make  them  applicable  to  multidimensional  cross- 
classifications. In  particular,  generalizations  of  the  linear-by- 
1 inear  association  model,  the  additive  effects  model  (model  2.6)  and 
the  multiplicative  effects  model  (model  2.5)  for  three  dimensional 
cross-classifications  of  ordinal  variables  are  discussed  in  the  first 
two  sections.  If  the  ordering  of  the  categorical  variables  is  taken 
into  account,  then  it  is  possible  to  construct  a richer  variety  of 
parsimonious  models  than  possible  with  the  standard  logl inear  model 
that  allows  partial  association  between  each  pair  of  variables.  Further- 
more, we  are  able  to  propose  unsaturated  models  that  allow  three-factor 
interaction  among  the  variables.  In  section  3 we  describe  briefly  how 
the  ML  estimates  of  the  model  parameters  can  be  obtained  using  the 
techniques  of  the  previous  chapter.  An  example  illustrating  the  fitting 
of  these  models  is  then  given  in  section  4.  Generalizations  for  com- 
binations of  nominal  and  ordinal  variables,  and  extensions  for  several 
dimensions  are  discussed  in  sections  5 and  6 respectively.  Finally,  we 
comment  on  model  selection. 

Much  of  the  material  presented  here  parallels  that  of  the  recent 
work  on  "partial  association"  models  of  Clogg  (1982)  and  Goodman  (1981). 
Whereas  they  define  models  via  properties  of  relevant  odds  ratios,  we 
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rather  stress  the  model  forms  for  log  expected  frequencies  in  an  attempt 
to  clarify  how  these  models  are  related  to  the  standard  log! inear 
models.  In  addition,  Clogg's  partial  association  models  do  not  permit 
three-factor  interaction. 

Consider  a three-way  table  with  r rows,  c columns  and  i layers. 

Let  X,  Y,  and  Z denote  the  row,  column  and  layer  variables  respectively. 

Let  m^k  denote  the  expected  frequency  under  the  stated  model, at  level 
i of  X,  level  j of  Y,  and  level  k of  Z,  for  l<i<r,  l<j<c,  and  l<k<£. 

To  describe  conditional  association  between  a pair  of  variables,  say 
row  and  column,  within  a fixed  level  of  the  third  variable  Z,  say  level 
k,  we  use  the  following  type  of  odds  ratios  referred  to  as  local  odds  ratios 

6ij(k)  = n,i jkmi+1 , j+1  ,k/mi+l . j,kmT,j+l  ,k  ls1sr-1-  1SJSC-1-  (3-D 

The  variables  X and  Y are  conditionally  independent  within  the  levels 
of  Z if  all  £(r-l ) (c-1 ) odds  ratios  equal  one.  If  the  level  of  the 
column  variable  Y or  the  row  variable  X is  fixed,  similar  odds  ratios 
9i(j)k  and  e(i)jk  describe  the  conditional  association  for  the  other 
pairs  of  variables.  Ratios  of  odds  ratios  {0^^}  are  used  to  describe 
the  local  three-factor  interaction,  where 

6ijk  " 6 i j ( k+1 )/0i j(k)  = 9i (j+1 )k/ei (j)k  = 6(i+l )jk/e(i )jk  * (3*2) 

There  is  no  three-factor  interaction  between  the  variables  if  all 
(r-1 ) (c-1 ) (£-1 ) of  the  0.jk  equal  one. 

3.1.  Additive  Effects  For  Three  Dimensions 

In  this  section  we  assume  that  the  table  has  all  categorical  variables 
ordered,  and  that  we  choose  (a  priori)  the  scores  {u. ,i=l,...,r}. 
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(Vj , j=l , . . ,c},  and  {wk,k=l , . . ,£}  to  assign  respectively  to  the  levels 
of  the  categorical  variables  X,  Y,  and  Z. 

3.1.1.  Homogeneous  Linear  Effects 

The  mutual  independence  model,  denoted  by  (X,Y,Z),  expresses  the 
expected  frequencies  {m..^}  as 

1 °9  mijk  = y + XX  + XY  + XZ,  where  ZXX  = IX  Y = IXk  = 0 . 

A parsimonious  logl inear  model  that  utilizes  the  ordinal  nature  of  the 
variables  and  allows  association  between  each  pair  of  variables  can  be 
expressed  as 

log  mij-k  = y + XX  + Xj  + Xk  + 3XY (u^u) (v^-V)  + BXZ(u. -u) (wk~w) 

+ BYZ(v-v)(wk-w)  . (3.3) 

This  model  is  a simple  generalization  of  the  independence  model,  with 
only  three  more  parameters.  It  has  residual  df  = rc£  - r - c - a - 1 , 
and  it  is  always  unsaturated. 

From  (3.3)  we  see  that 

109  eiJ(kf  tui+Tu1)(vj+Tvj)e3<Y 

108  ei(j)k = <“1+,->'1)<vrwk)exz 

109  e(i)jk = lVrV(Vi'“tl*,z 

log  61jk  = 0 . (3.4) 

Model  (3.3)  is  an  extension  of  the  1 inear-by-1 inear  association  model 
(2.4)  designed  for  two-way  tables.  In  other  words,  for  a fixed  level 
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of  a certain  variable  the  1 inear-by-1 inear  association  model  fits  the 
corresponding  two-way  cross-classification  formed  by  the  other  variables. 

In  addition,  the  same  association  parameter  applies  to  each  level  of  the 
third  variable.  The  model  assumes  no  three  factor  interaction,  since 
all  the  local  associations  are  homogeneous  across  the  levels  of  the  third 
variable.  Model  (3.3)  is  referred  to  as  a homogeneous  linear  effects  model. 

For  the  special  case  where  the  scores  are  equally  spaced,  the  local 
associations  described  in  (3.4)  are  constant  for  each  partial  table. 

For  instance,  for  integer  scores  {ui=i},  {v^j},  and  {wk=k},  the  local 
odds  ratios  satisfy  log  0.j(k)  = BX\  log  ei(j)|<  = BXZ,  and  log  e(.)jk  = BYZ. 
We  then  refer  to  the  model  as  a homogeneous  uniform  association  model. 

For  instance,  for  each  fixed  level  of  the  layer  variable  Z,  the  local 
odds  ratio  is  uniform  for  the  row  variable  X and  column  variable  Y,  and 
the  association  is  homogeneous  for  all  the  levels  of  the  variable  Z. 

The  values  of  the  association  parameters  are  dependent  on  the  choices 
of  categories  and  the  scores  assigned  to  those  categories.  To  reduce 
this  dependence,  we  can  rescale  the  chosen  scores  to  have  unit  standard 
deviation  in  each  dimension.  Then,  for  instance,  BXY  refers  to  con- 
ditional log  odds  ratios  for  levels  of  X that  are  one  standard  deviation 
apart  and  levels  of  Y that  are  one  standard  deviation  apart. 

If  model  (3.3)  fits  reasonably  well,  it  is  very  easy  to  describe 
the  pattern  of  association  in  the  table.  If  a particular  B»  say  BXY, 
equals  zero,  then  there  is  conditional  independence  between  X and  Y 
within  the  levels  of  Z.  Thus,  given  that  model  (3.3)  fits,  a single 
degree  of  freedom  test  of  conditional  independence  can  be  based  on  the 
difference  between  the  likelihood-ratio  chi-square  (G2)  for  model  (3.3) 
with  B =0  and  for  model  (3.3)  with  unrestricted  BXY  value 
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3-1 .2.  Homogeneous  Additive  Effects 

If  at  least  one  of  the  bivariate  conditional  associations  is  non- 
homogeneous,  other  effects  may  be  introduced  into  the  model  without 
affecting  the  structure  of  no  three  factor  interaction.  An  extension 
of  the  R+C  model  (2.6)  is  given  by 


log  m.jk  = y + + (uru)(Vj.-v)BXY  + (u.-uKw^w^ 


-N  „XZ 


— \„YZ 


+ (vj-v)  (w^-wJb14,  + + (wk-w)BXi  + (iij-u)B^j 

+ (wk-^)62j  + (ui_ir)6U  + (vj"V)62k’  (3.5) 


where  ZAv1 


nJ  - nl  * 0 a"d  6fl  ■ B?r  - ®2i  - ^2r  ' 8ll  * B,Yc  - bJ, 


Y Z Z Z Z 
" 62c  = 611  ‘ BU  " 621  " Bo.  = 0 . 


J2l 


This  model  has  an  additional  2[(r-2)  + (c-2)  + (£-2)]  linearly  independent 
parameters  besides  those  in  the  homogeneous  linear  effects  model  (3.3), 
so  df  = rc£  - 3r  - 3c  - 3£  + 1 1 . As  for  the  case  of  two-way  tables,  the 
constraint  on  the  B-parameters  is  clarified  when  we  look  at  the  local 
log  odds  ratios. 

From  (3.5)  we  see  that 


109  eij(k)  - (Vrui)<vj+rVB  + (Vrvj)(Bi,i+r6i,i> 


log  eijk  = o . (3.6) 

For  the  special  case  of  equal -interval  scores,  the  local  log  odds  ratio 
takes  the  form 

109  6ij(k)  = 6 + + @5  • 


(3.7) 


35 


We  remark  from  (3.6)  and  (3.7)  that  the  local  odds-ratios  {0. 

i j vkj 

have  the  same  form  as  the  R+C  model  (2.6)  does.  That  is,  for  each 
fixed  level  of  the  layer  variable  Z,  the  local  log  odds-ratios  are 
linear  functions  of  the  differences  in  the  relevant  row-effects  and 
column-effects.  Also,  these  additive  effects  are  homogeneous  across 
the  levels  of  Z.  Hence,  we  refer  to  model  (3.5)  as  a homogeneous 
additive  effects  model.  When  all  6^.  = 0,  we  then  obtain  a model  with 
homogeneous  additive  column  effects  (i.e.  log  e.^j  = a.jbjBXY  + c^Bj) 
on  the  X-Y  association.  In  other  words,  the  same  C model  (2.3)  is  fitted 
to  each  level  of  the  variable  Z.  When  all  B^-  = 0 and  BYj  = 0,  for 
integer  scores  we  obtain  homogeneous  uniform  X-Y  association  e..,  . = 

Y V 1 J \ "W 

exp(B  )•  Generally,  the  pattern  in  the  (BX  and  (b|  j^-B^} 

indicate  how  the  conditional  X-Y  local  odds  ratios  depart  from 
uniformity. 

3.1.3.  Interaction  Models 

If  models  (3.3)  and  (3.5)  fail  to  fit,  either  there  is  three- 
factor  interaction  or  else  the  homogeneous  association  effects  are 
not  the  simple  additive  ones  of  model  (3.5).  In  this  subsection  we 
consider  generalizations  of  models  (3.3)  and  (3.5)  that  allow  for 
three-factor  interaction. 

A simple  generalization  of  the  homogeneous  linear  effects  model 
(3.3)  that  allows  for  three-factor  interaction  is  obtained  by  adding  a 
single  parameter  for  interaction.  The  model 

109  mijk  = V>  + ^ + + xk  + (Vu)(vj-V)BXY  + (ui-u)(w|<-w)BXZ 

+ (Vj-v)(wk-w)BYZ  + (u^-IT)  (Vj-7)  (wk~w)BXYZ 


(3.8) 
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has  residual  df  = rcl  - r - c - i - 2.  For  this  model. 


-v,)[SXY  ♦ (wk-w)6XVZ] 


(3.9) 


From  (3.8)  and  (3.9)  we  see  that  (i)  the  departure  of  log  m. ..  from 

1 J K 


independence  is  a linear  function  of  each  categorical  variable  given  that 
the  other  variables  are  fixed  at  each  level,  and  (ii)  the  association 
between  any  two  variables  is  nonhomogeneous  across  the  levels  of  the 
third  variable  but  changes  linearly.  The  model  is  referred  to  as  the 
linear  association  effects,  linear  interaction  effects  model.  The  inter- 
action is  constant  for  all  2x2x2  subtables  with  the  same  distance 
between  rows,  the  same  distance  between  columns,  and  the  same  distance  between 
layers.  For  integer  scores,  the  local  interaction  is  constant,  and 
model  (3.8)  can  therefore  be  described  as  a uniform  interaction  model. 

In  that  case,  we  note  that  from  (3.9)  the  association  between  X and  Y 
is  uniform  within  each  layer  of  Z,  but  the  strength  of  association  is 
a linear  function  of  Z.  Hence  model  (3.8)  is  a type  of  heterogeneous 
uniform  association  model . 

When  the  association  between  each  pair  of  variables  is  not  homo- 
geneous within  the  levels  of  the  third  variable,  we  can  also  construct 
a model  that  allows  additive  effects  but  maintains  constant  local  inter- 
actions. This  model  can  be  obtained  by  combining  model  (3.3)  and  model 
(3.8).  Specifically,  we  add  the  interaction  term  B^(u .-17) (v  .-V) (w.  -w) 
to  model  (3.3).  Hence  the  resulting  model  contains  row,  column,  and 
layer  effects  on  the  partial  associations,  and  again  linear  heterogeneity 
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in  the  level  of  the  local  conditional  associations.  That  is,  for 
instance,  the  X-Y  association  changes  linearly  across  the  levels  of  Z. 
We  refer  to  the  model  as  a homogeneous,  additive  association  effects, 
linear  interaction  effects  model . 

If  any  of  the  row,  column,  or  layer  effects  are  not  homogeneous, 
the  interaction  will  not  be  uniform.  A general  model  of  this  type, 
which  utilizes  the  ordinal  nature  of  the  variables,  is 


1og  mijk  = y + Ai  + xj  + Ak  + (ui-u)(Vj-v)(3XY+8*) 

+ (ui-u)(wk-w)(6XZ+6j)  + (vj.-v)(wk-w)(BYZ+6-) 

+ + (w^w^  + ( u i -T7)  + (wk~w)  62  j 

+ (Ui-u)BZk  + (Vj-V) 3|k  , (3.10) 

where  BX  = 3*  = $Y  = BY  = BZ  = BZ  = 0 . 

For  this  model,  df  = rc£  - 4(r+c+£)  + 16.  From  (3.10)  we  obtain  for  the 
case  of  integer  scores 


log  6 


log  9 


ij(k)  = exv  ^*(®{,i+r»!i)Me;,i+1-8;j) 

+ (k-(J.+l)/2)[BXYZ+(6X+1-6X)  + (6Y+1-6Y)]  , 
m - + (Bi+r6i> + (ej+r6j> + • 


(3.11) 


Note  that  the  row,  column,  and  layer  effects  change  linearly  over  the 
levels  of  other  variables.  The  interaction  is  a linear  function  of  the 
differences  between  row,  column,  and  layer  effects.  A variety  of 
structures  are  special  cases  of  this  model,  including  all  models  con- 
sidered previously.  We  refer  to  the  model  as  the  additive  interaction  model. 
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In  some  cases,  we  may  wish  to  let  some  of  the  effects  be  as  general 
as  for  nominal  variables.  In  such  a case,  depending  on  the  aim  of  the 
researcher,  more  general  models  can  be  constructed.  For  example, 
suppose  that  our  main  interest  is  in  the  X-Y  conditional  association, 
and  that  we  posit  heterogeneous  linear-by-1 inear  X-Y  association  but 
have  no  preconceived  notions  about  the  natures  of  the  Y-Z  and  X-Z 
associations.  Then,  we  could  fit  the  model 


XZ  A ,YZ 


109  mijk  = y + Xi  + Xj  + Ak  + Xik  + Xjk  + (Vu)(Vj-v)(eXY+BkYZ)  (3.12) 


XYZ 

with  £Bk  = 0,  which  has  residual  df  = £(rc-r-c).  Fitting  this  model 
corresponds  to  fitting  the  bivariate  1 inear-by-1 inear  association  model 
(2.4)  separately  at  each  level  of  Z.  One  set  of  likelihood  equations 
satisfied  by  this  model  is  that  the  observed  X-Z  marginal  counts  equal 
the  fitted  X-Zmarginal  counts  (i.e.  fr,i+k  = ni+k  l^i<r).  Hence,  we  also 
might  use  a model  like  this  if  the  X-Z  or  Y-Z  marginal  table  is  fixed  by  the 
sampling  scheme.  It  is  the  general  heterogeneous,  1 inear-by-1 inear  X-Y 
association  model. 

The  most  complex  form  of  the  linear  interaction  model  has  the  form 


XY 


. XZ  . , YZ 


log  mij(<  = v + \*  + \ . + + A,,  + + A, 


ij 


Ik  + xjk  + <Vu)(vj-v)(Vw)B 


— \nXYZ 


(3.13) 


This  is  the  general  (XY,XZ,YZ)  loglinear  model,  with  structured  inter- 
action term  added.  This  model  has  residual  df  = (r-1 ) (c-1 ) (£-1 )-l . 
When  applied  with  integer  scores,  it  is  the  uniform  interaction  model 
suggested  by  Goodman  (1979,  1981). 

Given  that  a model  such  as  (3.8)  or  (3.13)  fits,  a single  degree- 
of-freedom  chi-square  test  of  no  three-factor  interaction  can  be  based 
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2 

on  the  difference  between  the  G values  for  the  model  with  and  without 

the  interaction  term.  Given  that  a model  such  as  (3.12)  holds,  a 

df  = SL  - 1 test  of  independence  of  no  interaction  can  be  based  on  the 

2 

difference  between  the  G value  for  the  model  with  heterogeneous 

XY  XY7 

association  parameters  (p  + Bk  L)  and  the  model  with  homogeneous 

XY 

association  parameter  B • 

Notice  that  we  have  formulated  these  models  according  to  the 

hierarchy  principle  followed  by  standard  loglinear  models.  Presence 

of  higher  order  interaction  terms  implies  that  the  lower  order  terms 

should  be  present  having  the  same  or  more  complex  effects.  For  example, 

it  would  usually  not  be  of  interest  to  have  a general  interaction  term 
XYZ 

^ijk  ^near  association  effects.  Figure  1 illustrates  the  relation- 
ships among  the  models  we  have  discussed  in  this  section. 

3.2.  Multiplicative  Effects  For  Three  Dimensions 

In  the  RC  model  (2.5)  the  scores  for  the  levels  of  the  ordinal 
variables  are  estimated  rather  than  fixed.  An  equivalent  model  to  the 
homogeneous  linear-effects  model  (3.3)  can  be  expressed  as 

l°g  m1jk  = p + AX  + AY  + AZ  + e^Vj  + 8XZUjMk  + SYZvjB|c  , (3.14) 

where  IAX  = EAY  = £Ak  = Zpf  = = Zn)k  = 0 and  Eu?  = Zv?  = = 1 . 


We  have  taken  the  basic  independence  model  and  added  three  association 
parameters  and  (r-2)  + (c-2)  + (£-2)  independent  score  parameters,  so 
df  = rc£  - 2(r+c+£)  + 5.  This  is  always  an  unsaturated  model,  and  it 
is  easy  to  interpret  since 


,XY 


log  ejj(k)  = b (gi+1-v1)(vj+1-vJ)  , 


'og  e1jk 


= 0 . 


(3.15) 
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Figure  1.  Additive  effects  model  for  ordinal  variables  (arrows 
indicate  generalizations  of  models). 
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In  fitting  this  model  we  are  estimating  scores  under  which  the 
homogeneous  linear  effects  property  holds.  That  is,  distances  between 
estimated  scores  indicate  what  the  spacings  between  rows,  between 
columns,  and  between  layers  must  be  so  that  the  partial  odds  ratios 
for  two  variables  depend  only  on  distances  between  levels  and  is  the 
same  across  each  level  of  the  third  variable. 

Unlike  the  models  of  section  2,  model  (3.14)  is  not  loglinear  in 
its  natural  parameters,  and  it  is  invariant  to  re-orderings  of  the 
categories.  Hence,  the  parameter  scores  need  not  have  the  same 
ordering  as  the  levels  of  the  ordinal  variables.  If  the  model  fits, 
lack  of  monotonicity  in  the  scores  for  a variable  indicates  nonmonotonic 
associations,  in  the  sense  that  local  associations  involving  that 
variable  are  positive  in  some  locations  and  negative  in  others. 

A parameter-scores  version  of  the  linear  interaction  effects  model 
(3.8)  with  an  additional  parameter  6XYZ  beyond  these  in  model  (3.14),  is 

log  m1jk  = y + XX  + iY  + XZ  + „.VjeXY  + 'P^kBXZ  + vjukBVZ 

XY7 

+ ^ivj^ke  (3.16) 

Here,  df  = rci,  - 2(r+c+£)  + 4,  and 

109  6ij(k)  = (yi+l'u1  )('Jj+r%Jj)[6XY  + (aik-u)BXYZ] 

log  eijk  = (y,+1-P1)(Vj+1-vJ)(uk+1-iiik)6XYZ  ..  (3.17) 

Treating  the  },  {vj},  and  {a^}  as  scores,  we  see  that  we  can  interpret 
model  (3.16)  as  a heterogeneous  linear  association  effects  model  that 


has  a linear  interaction. 
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The  estimated  scores  might  depart  very  much  from  monotonicity. 

In  that  case,  it  makes  more  sense  to  interpret  them  as  row,  column,  and 
layer  effects.  Model  (3.16)  can  be  viewed  similarly,  but  with  hetero- 
geneity in  the  overall  levels  of  bivariate  association. 

Models  (3.14)  and  (3.16)  can  be  further  generalized  in  obvious 
ways.  For  example,  greater  flexibility  is  provided  by  letting  the 
effects  for  each  variable  differ  for  each  association,  as  in  the  model 


log  m 


i jk 


= y + X?  + 


+ X 


+ ylivlj 


,XY 


y2iV 


,6 


XZ 


v2jW2kB 


YZ 


(3.18) 


This  model,  which  has  df  = rc£  - 3(r+c+£)  + 11,  is  the  multiplicative- 
effects  version  of  model  (3.5).  We  refer  to  it  as  the  homogeneous 
multiplicative  effects  model . 

A natural  interaction  term  for  model  (3.18)  has  the  multiplicative 
Xyz 

form  as  in  model  (3.16).  With  this  generalization,  this 

model  can  be  interpreted  as  one  having  homogeneous  multiplicative  row, 
column,  and  layer  effects,  but  heterogeneity  in  the  overall  levels  of 
conditional  association  since  log  has  the  form  ei  e . e 

The  model  is  referred  to  as  the  heterogeneous  multiplicative  effects 
model . 


3.3.  Maximum  Likelihood  Estimation  of  Model  Parameters 

The  estimates  of  the  model  parameters  can  be  obtained  through  the 
use  of  the  existing  available  computing  procedures  described  in  Chapter 
Two.  To  illustrate,  let  us  consider  the  additive  effects  model  (3.3) 
with  the  addition  of  a linear  interaction  term,  i.e.,  the  composite  of 

models  (3.3)  and  (3.9).  The  kernel  of  the  log  likelihood  in  that  case 
is 
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. En.ijk  l09  mijk  = N|J  + &W?  + £n+J+*j  + £Wk 

1 j J » K 


,XY 


+ BA'  E n.j+(u.-u)(vj.-v)  + 3XZ  Zni+k(u.-u)(wk-w) 

i>J  i,k 


+ 6XYZZ(u1-u)(vj-v)(wk-w)niJk 


(3.19) 


Taking  the  partial  derivatives  of  the  loglikelihood  with  respect 
to  all  the  parameters  and  setting  them  equal  to  zero,  we  obtain  the 
following  set  of  equations  to  solve: 


A 


mi++ 


n 


i++ 


A 


m++k  “ n++k 


3 Aj+(xru)(vk-v)  = 

* 9 J 


S nij+(u.-u)(vrv) 

* 9 J 

/kni+k(ur“><Vw> 


(3.20a) 

(3.20b) 

(3.20c) 

(3.21a) 


(3.21b) 
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ASi+jk(vrv)(wk't')  * ? tn+jk(vj-v)(wk-w) 

J » ^ J 9 K 

JSjj+(uru)  = zn1j+(uru) 


(3.21c) 


(3.22a) 


^i+k(ui‘u)  = f i+k(ui-u> 


(3.22b) 


fij+(vj'V)  = fij+<Yv) 

J J 


(3.22c) 


?3+ok(vrv> = f+jk<vj-v) 


(3.22d) 


jfl+k(Vw>  * j;Vk(Vw> 


(3.22e) 


^jk<Vw>  ■ sn+3k(wk-w) 


(3.22f) 


£ i«ijk(ui-u)(vj-v)(vw> = . 5 "ijk(uru)<vj-v)K 

1 I * J » K 


-w)  (3.23) 


Equations  (3.21)  are  implied  by  (3.22),  so  they  are  redundant  for 
the  model  in  its  present  form.  Note  from  (3.23)  that  an  implication  of 
this  three-factor  interaction  model  is  that  the  observed  and  fitted  third 
moments  are  equal.  If  we  replaced  (w^-w)B^^  by  the  heterogeneity 
interaction  term  3^  , we  would  have  the  further  restriction  that  the 
observed  X-Y  correlation  equals  the  fitted  X-Y  correlation  within  every 
level  of  Z.  That  is,  ^.rn-j  j k (ui  “u ) (vj”v)  = 2 n ^ ^ k ( u ^ -TJ)  ( v j - V)  for  k = 

Subsets  of  the  above  constraints  are  satisfied  by  the  additive  models 
(3.3),  (3.5),  and  (3.8).  The  homogeneous  linear  effects  model  satisfies 
(3.20)  and  (3.21)  only.  That  is,  the  observed  and  fitted  frequencies 
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for  model  (3.3)  have  the  same  marginal  correlations.  This  is  also 
true  for  the  homogeneous  additive  effects  model  (3.5),  which  also 
satisfies  (3.22).  Equation  (3.22a)  implies  that  within  each  level  of 
Y,  the  observed  mean  of  X equals  the  fitted  mean  of  X (similar 
descriptions  apply  to  (3.22b )- (3- 22b) ) . The  linear  association-effects, 
linear  interaction-effects  model  (3.8)  satisfies  (3.20),  (3.21),  and 
(3.23).  That  is,  the  observed  and  fitted  marginal  correlations  are 
equal,  and  the  observed  and  fitted  third  moments  are  equal. 

As  in  the  previous  chapter,  the  logl inear  models  discussed  in 
sections  1,  5,  and  6 are  special  cases  of  Nelder  and  Wedderburn's 
(1972)  generalized  linear  models.  Hence  the  computer  package  GLIM  can 
be  used  to  fit  these  models.  The  asymptotic  distribution  of  the 
parameter  estimates  of  these  models  is  normal.  For  a proof  of  asymptotic 
normality  of  the  ML  estimates  of  parameters  of  models  designed  for  two- 
way  tables,  we  refer  the  reader  to  Appendix  1.  A generalization  of 
that  is  then  obvious. 

The  computer  package  GLIM  produces  the  ML  estimates  of  the  model 
parameters  and  the  asymptotic  covariance  matrix  of  the  estimates. 
Alternatively,  the  LOGLINEAR  procedure  contained  in  the  computer 
package  SPSSX  can  be  used  for  this  purpose.  Although  the  two  existing 
procedures  above  are  designed  for  logl inear  models,  they  still  can  be 
used  to  find  the  ML  estimates  of  the  parameters  of  multiplicative 
effects  models.  This  can  be  done  by  estimating  only  one  set  of  parameter- 
scores  at  a time,  and  treating  the  remaining  sets  as  fixed.  Then  the 
process  is  as  described  in  the  previous  chapter. 
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3.4.  Example 

The  data  in  Table  3.1,  taken  from  the  book  by  Fienberg  (1980,  p.  6), 
were  originally  reported  by  Cornfield  (1962).  The  sample  consisted  of 
1329  male  residents  of  Framingham,  Massachusetts,  aged  40-59,  classi- 
fied by  blood  pressure  (P)  and  serum  cholesterol  (C)  levels.  During  a 
six-year  follow-up  period,  they  were  also  classified  according  to 
whether  they  developed  coronary  heart  disease  (H).  The  categories  of 
each  variable  are  ordered,  and  we  will  use  Table  3.1  to  illustrate  the 
ordinal  models  discussed  previously. 

These  data  could  be  analyzed  in  several  ways.  For  instance. 

Cornfield  applies  a logit  model  which  treats  H as  a response  variable. 
However,  for  the  association  models  in  this  chapter,  we  need  not 
identify  a response  variable,  and  we  also  obtain  information  about  the 
C-P  association.  Table  3.2  contains  the  results  of  fitting  four  of 
these  models  to  Table  3.1: 

(a)  Model  (3.3)  with  integer  scores  (homogeneous  uniform  association) 

(b)  Model  (3.5)  with  integer  scores  (homogeneous  additive  effects) 

(c)  Model  (3.14)  (parameter-scores  version  of  (3.3)) 

(d)  Model  (3.18)  (homogeneous  multiplicative  effects) 

Models  (3.5)  and  (3.14)  are  generalizations  of  (3.3),  and  model  (3.18) 

is  a generalization  of  (3.14).  All  four  models  fit  adequately.  The 
2 

decrease  in  G obtained  with  the  more  complex  models  is  only  on  the 
order  of  the  decrease  in  degrees  of  freedom.  In  order  to  select  a 
model  for  Table  3.1,  we  must  balance  the  parsimony  of  structure  and 
interpretation  of  a model  versus  the  complexity,  but  more  complete 
sample  descriptions  of  another  model.  For  instance,  it  is  easy  to 
interpret  model  (3.3),  whereas  models  (3.5)  and  (3.18)  are  very  complex, 
but  provide  more  complete  descriptions  of  the  data. 


47 


Table  3.1.  Cross-classification  of  Framingham  men  by  blood  pressure, 

serum  cholesterol,  and  presence  or  absence  of  heart  disease. 


Coronary 

Serum 

Systolic 

Heart 

Cholesterol 

Blood  Pressure  (mm  Hq) 

Disease 

(mg/100  cc) 

<127 

127-146 

147-166 

167+ 

<200 

2 

3 

3 

4 

200-219 

3 

2 

0 

3 

Present 

220-259 

8 

11 

6 

6 

>260 

7 

12 

11 

11 

<200 

117 

121 

47 

22 

200-219 

85 

98 

43 

20 

Absent 

220-259 

119 

209 

68 

43 

£260 

67 

99 

46 

33 

Source: 

Fienberg  (1980) 
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According  to  the  simple  uniform  association  model,  the  data  may 
be  well  described  by  constant  conditional  log  odds  ratios  of  .53  for 
the  H-C  association,  .44  for  the  H-P  association,  and  .10  for  the  C-P 
association.  All  three  association  terms  are  needed  in  the  model, 
since  their  respective  asymptotic  standard  errors  are  .11,  .12,  and 
.03.  Thus,  model  (3.1)  suggests,  for  instance,  whether  the  heart 
disease  is  present  or  not,  the  higher  the  level  of  serum  cholesterol 
in  those  men,  the  higher  their  blood  pressure  would  tend  to  be.  Also, 
the  positive  C-P  association  is  uniform  across  the  levels  of  coronary 
heart  disease.  Similar  descriptions  apply  for  any  pair  of  variables 
across  the  levels  of  the  third  variable.  To  make  the  association 
parameter  estimates  for  the  uniform  association  model  comparable  in 
value  to  those  in  the  multiplicative  models  of  similar  form,  we  can 
rescale  the  fixed  scores  to  satisfy  the  constraints  Eu?  = Zv2  = Ew2  = 1 . 

*1  J K 

The  association  parameter  estimates  for  model  (3.3)  then  equal 
~HC  ~HP  -■'CP 

8 - «84,  8 - .70,  and  8 = .48,  very  similar  to  those  of  models 

(3.14)  and  (3.18). 


For  illustrative  purposes  we  will  discuss  the  other  models,  even 
though  they  do  not  suggest  significant  departures  from  uniform 
associations.  Let  (i,j,k)  denote  level  i of  H,  level  j of  C,  and 
level  k of  P.  For  the  homogeneous  additive  effects  model  (3.5), 
differences  in  8-parameters  measure  departures  of  conditional  pairwise 
association  from  uniformity,  as  described  from  (3.7).  Since  the  local 
log  odds  ratio  log  0^.^  = §HP  + (8^-8^)  = -.07  is  close  to  zero 
compared  to  log  © i ( j ) 2 anc*  ^l(j)3’  the  layer  effect  estimates 

^ P 

Pk^  indicate  that  the  first  two  levels  of  P have  similar  distributions 
with  respect  to  H in  the  sample.  Similarly,  the  column  effect 
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Table  3.2.  Results  of  fitting  association  models  to  Table  3.1. 


Goodness- 

of-fit 

(3.3) 

Homogeneous 

Uniform 

Association 

(3.5)  Mode1  (3.14) 
Homogeneous  Homo.  Linear 
Additive  Effects 

Effects  Para.  Scores 

(3.18) 
Homogeneous 
Mul tipi icative 
Effects 

G2 

22.8 

13.6 

18.5 

13.5 

df 

21 

13 

17 

13 

Parameter  Estimates 

H-C 

SHC=  .53 

BHC=  .45 

BHC=  .82 

6HC=  .86 

Association 

8^2=  -.66 

y-j  - -.46 

-C  _ ^ r 

Ul  1 - -.35 

§i3=  -.36 

~c_  . « 

y2-  -.42 
r 

.C  _ r-| 

U12-  *.51 

r 

\if  .10  y^=  .08 

P$=  .78  yf4=  .78 


H-P 

§HP- 

.44 

eHP= 

.40 

8HP= 

.70 

bhp= 

.73 

Association 

gP  - 
B12 

-.47 

~P 

V1  - 

-.60 

~p  . 

vir 

-.41 

*P  = 

B13 

-.27 

~P 

v2  = 

-.27 

.p 

V1 2 

-.48 

~P 

v3  " 

.13 

~P 

v13" 

.12 

- 

v4  " 

.74 

-P  . 
v14" 

.77 

C-P 

gCP= 

.10 

BCP= 

.10 

BCP= 

.48 

Bcp= 

.46 

Association 

B2C2= 

*23= 

§P  = 
S22 

6P  = 
B23 

-.03 

-.04 

.07 

-.01 

-c  _ 

yl  - 

y2  ' 
-C  _ 
y3  " 

y4  " 

^ D 

-.46 

-.42 

.10 

.78 

y21~ 

-C 

y22 

y23 

y24" 

-.62 

-.34 

.37 

.60 

Ar 

vi  = 

D 

-.60 

~P 

v21~ 

-.78 

/V  1 

v2  = 
P 

-.27 

-P 

v22 

.11 

A 1 

v3  = 

^ P 

.13 

„P 

v23 

.06 

^ r 

v4  = 

.74 

~P 

v24 

.61 

Note:  For  homogeneous  additive  effects  model,  =^i4=^i i =BP4=0, 
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Table  3.3.  Estimated  local 

log  odds  ratios 

for  models 

fit  to  Table  3.1 

Log  odds 

(3.3) 

Homogeneous 

Uniform 

(3.5)  Model  (3.14) 
Homogeneous  Homo.  Linear 
Additive  Effects 

(3.18) 

Homogeneous 

Multiplicative 

ratio 

Association 

Effects  Para.  Scores 

Effects 

H-C  Association: 

i°g  §11(k) 

.53 

-.20 

.05 

-.19 

log  e12(k) 

.53 

.75 

.60 

.73 

1og  ®13(k) 

.53 

.81 

.79 

.85 

H-P  Association: 

l09§i(j)i 

.44 

-.07 

.32 

-.07 

log  01(;j)2 

.44 

.60 

.40 

.62 

los  @l(j)3 

.44 

.72 

.61 

.67 

C-P  Association: 

log  §(1))1 

.10 

.14 

.01 

.11 

log  e(j)2, 

.10 

.16 

.08 

.29 

103  §(i)31 

.10 

.21 

.11 

.09 

109  ®(i)12 

.10 

-.01 

.01 

-.01 

log  e(j)22 

.10 

.01 

.10 

-.01 

109  §(i)32 

.10 

.06 

.13 

-.01 

109  ®(i)13 

.10 

.08 

.01 

.07 

109  S(i)23 

.10 

.10 

.15 

.18 

109  9(i )33 

.10 

.15 

.20 

.06 
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^ P 

parameters  (B-j'j}  indicate  that  the  first  two  levels  of  C have  similar 
distributions  on  H,  with  log  8^^  = -.20  representing  a slight 
tendency  in  the  sample  for  the  incidence  of  heart  disease  to  be  lower 
at  the  second  cholesterol  level  than  at  the  first.  The  and 

1 S2 j ) indicate  that  the  middle  two  levels  of  P have  similar  distributions 
with  respect  to  C. 

For  the  case  when  the  scores  are  estimated  rather  than  selected  in 
advance,  we  reach  the  same  conclusions  as  those  by  previous  models.  For 
the  homogeneous  linear  effect  parameter  scores  model  (3.14),  the  score 
estimates  are  monotonic  but  ylj'  and  are  quite  close,  indicating  that 
the  first  two  levels  of  C had  similar  distribution  on  H,  and  similar 
distributions  on  P in  the  sample.  For  the  more  complex  model  (3.18), 
we  permit  the  effects  (or  scores)  for  each  variable  differ  for 
each  partial  association.  For  instance,  if  we  consider  the  H-C 
association  the  estimated  parameters  y^  and  y^  are  still  quite  close 
indicating  that  the  first  two  levels  of  C had  similar  distributions  on 
H,  but  since  y-j^y-j^,  it  makes  more  sense  to  interpret  them  as 
cholesterol -effects  than  scores.  Similar  remarks  can  be  drawn  for  the 
H-P  association.  Departures  from  uniform  association  are  the  same  as 
the  ones  suggested  by  model  (3.5).  The  sample  data  suggest  through 
both  models  (3.5)  and  (3.18)  that  there  may  be  some  threshold  below 
which  C and  P are  not  associated  with  H.  For  further  details  on  the 
association  patterns  suggested  by  the  four  ordinal  models,  see  Table 
4.3,  which  gives  the  estimated  local  log  odds  ratios  for  each  partial 


association. 
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3.5.  Generalizations  For  Several  Dimensions 

Suppose  now  that  we  have  a cross-classification  of  d ordinal 

variables  (X^  ,X2, . . . ,Xd),  with  preassigned  scores  (uj1-,  u^) 

at  the  i.  levels  of  the  ith  variable,  i = 1,2 d.  Let  m.  denote 

the  expected  frequency  in  cell  j_=  ( i -j , i 2 , — * i d ) • The  independence 
model  is 

d X 

log  m.  = y + E x.d  . 

1 j=i  j 

A simple  extension  of  this  model  which  utilizes  the  ordinal  nature  of 
the  variables  is 

X XX 

log  m.  = y + IX. J + z (u!^-u^^)(u(.k^k^)6  J k (3.24) 

J J j<k  \ 


This  model  is  an  extension  of  the  homogeneous  linear  effects  model  for 

the  case  in  which  the  cross-classification  has  arbitrary  dimensions.  It 

states  that  there  is  a homogeneous  1 inear-by-1 inear  conditional  association 

X.X. 

between  each  pair  of  variables,  with  6 J being  the  log  odds  per  unit 

distances  on  X^.  and  X^  within  all  combinations  of  fixed  levels  of  the 

other  variables.  Model  (3.24)  has  an  additional  d(d-l)/2  parameters 

d 

besides  those  in  the  independence  model,  so  df  = - [1+  z (£.-l) 

12  d ._.| ' i ' 

+ d(d-l  )/2] . If  we  treat  the  variables  as  nominal  and  fit  the  standard 

logl inear  model  {X-jX^jX-jX^, . . . , Xcj_i X^ } having  all  bivariate  partial 

d 

associations,  then  the  residual  df  = [1  + z U--1)  + Z (l  .-1 ) (z.  -1 )] 

i=l  1 j>k  J k 


can  be  considerably  less  than  the  residual  df  of  the  proposed  model, 
especially  when  a couple  of  the  variables  have  several  categories. 
Model  (3.24)  can  be  generalized  by  adding  some  or  all  the  (!?) 

O 

additional  three  factor  interaction  terms  z (u(^-u^J ^ ) (u ^-IT^ ) 

j<k<£  nj  ik 
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/ (j?,)  — (o). 

(u.  -uv  ')S  J . For  this  three  factor  interaction  model,  the 

A 

association  between  each  pair  of  variables  is  homogeneous  within  levels 
of  other  variables,  but  each  log  odds  ratio  changes  linearly  across 
the  levels  of  each  of  the  other  variables. 

A parameter-scores  version  of  model  (3.24)  expresses  the  expected 
frequencies  as 


log  m.  = y + + Z yPV^e  J k . 

1 J ’j  j<k  nj  \ 


(3.25) 


This  model  provides  greater  flexibility  by  allowing  for  the  estimation 
of  scores  under  which  the  homogeneous  linear  effects  model  fits  best. 

The  model  has  residual  df  = A^...^  - [1  + z(A.-l)  + z(A.-2)  + d(d+l)/2]. 

Model  (3.24)  can  be  generalized  by  adding  homogeneous  or  heterogeneous 
variable  effects  (as  in  model  (3.10)  for  the  case  d = 3).  However,  the 
number  of  such  parameters  can  be  very  large  when  d>4,  and  the  resulting 
model  can  be  quite  complex  to  interpret. 


3- 6-  Generalizations  For  Combinations  of  Nominal  and  Ordinal  Variables 
Suppose  that  (X,,...,X.  ) are  ordinal  variables,  (X,  X.  ..  ) 

1 a]  d-j  + l ^l+°2 

are  nominal  variables,  and  let  d = d1+d2>  For  this  case,  a simple 
generalization  of  the  independence  model  which  allows  all  pairwise 
associations  but  no  three-factor  interaction  is 


log  m.  = 


y + + z A?k  + z (u!-i)-u(j))(u(k)-u(k))3XJXk 

J J d-j  < j<k  j'k  j<k<d1  j \ 

+ Z z (u(k)-u(k))3Xj(k)  , 
j>d1  k<d]  k 1j 


(3.26) 
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X-  XX  XX  X.(k) 

where  LA,J  = 0 all  j,  EA.J.  = ZA.J_.  = 0 all  j<k,  and  Z g.J  = 0 

ij  'o  ij  ’j’k  ik  ’j’k  ij  ’j 

x,xk 

for  all  j.  The  (£--l)(£^-l)  independent  parameters  (for  fixed 

J K ’j’k 


(j»k) ) pertain  to  the  associations  between  the  nominal  variables  X.  and 

Mk 

Xk’  di<J<k*  The  $ parameter  pertains  to  the  conditional  linear-by- 

1 inear  association  between  the  ordinal  variables  X.  and  X,  , j<k<d, , 

J K 

within  all  combinations  of  levels  of  the  other  variables.  The  (£.  ,) 

XJk)  J‘ 

independent  parameters  (for  fixed  (j,k))  represent  effects  of 

J 

the  levels  of  nominal  variable  X,  on  its  association  with  the  ordinal 

J 

variable  X^,  l^k^d^,  d ^ < j , like  the  association  term  in  the  R model 
for  the  case  d1  = 1 and  d2  = 1 . This  model  has 


df  = £,£«... £rf  - [1  + Z(£.-l)  + E (£,-l)(£.-l)  + (Jl ) 

j J d-j  < j<k  J k l 


+ d,  Z (£i-l)] 
1 d-|<j  J 


Model  (3.26)  tends  to  be  a much  simpler  and  more  parsimonious  model  than 

the  standard  logl inear  model  containing  all  bivariate  associations, 

especially  when  the  ordinal  variables  have  several  categories.  It  may 

XX 

be  easily  interpreted  through  the  odds  ratios  ead  K describing  local 
association  at  levels  a and  a+1  of  X.  and  levels  b and  b+1  of  X,  , within 

J K 

all  fixed  combinations  of  the  other  variables.  For  two  ordinal  variables, 
the  conditional  local  log  odds  ratios  equal 


X,X 


lo9 k ■ (uii,>-u<J))(u<|;).u(k))6xiXk , 


ab 


a+1  a 


b+1  ub 


Vk 


so  3 describes  the  local  log  odds  ratio  per  unit  distance.  For 
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integer  scores,  the  conditional  (Xj,Xk)  associations  are  characterized 
by  uniform  association.  For  X.  nominal  and  X.  ordinal,  we  have 

J K 

lo9  6ab  k = (ub+1-ubk))(Ba+l<k>-eaJ(k))  • 


'b+1  b 


a+1 


Thus,  for  equal  interval  scores,  the  difference  in  8-parameters 

describes  the  difference  between  the  two  levels  of  X.  with  respect  to 

the  ordinal  variable  on  the  log  scale.  Since  X.  is  nominal,  there  is 

J 

no  reason  for  selecting  adjacent  levels  of  X.  for  forming  these  odds 

J 

ratios. 

The  parameter-scores  analog  of  (3.26)  is  useful  for  providing 
more  flexibility  in  fitting  the  simple  type  of  nominal-ordinal  model. 

This  model  has  the  form 

X,  X.X.  X.X, 

log  m.  = y + ZA,J  + Z A • • + Z y>J  V*  B J * 

- J dj<j<k  j k j<k<d-|  \ 

+ Z Z y!k)6^'(k)  . (3.27) 

j>d]  k<d]  1k  nj 

The  model  has  Z fewer  degrees  of  freedom  than  does  model 

j<d1  J 

(3.26)  when  homogeneous  scores  are  used,  but  it  has  the  same  type  of 
simple  interpretation. 

To  illustrate  these  types  of  models  we  consider  the  case  d^  = 2 
and  d£  = 1 . That  is,  the  first  two  variables  X,Y  are  ordinal  whereas 
the  third  variable  Z is  nominal.  For  the  case  d = 3,  model  (3.26)  can 
be  expressed  as 


l°g  mij.k  - y + AC  + \\  + Ak  + (ui-u)(vj-v)8XY  + (u.-u)8(k 


-sJ. 


+ 


(Vj-vjB 


Z 

2k  ’ 
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where  ZAX  = ZXY  = ZAk  = Z6^k  = ZB^  = 0.  This  model  is  invariant  to 

re-orderings  to  the  levels  of  Z.  For  integer  scores,  the  model  implies 

uniform  conditional  X-Y  association  (log  e..,.,  = BXY),  log  odds- 
Z Z ^ 

ratios  of  B-|k.  - Blk  when  (within  levels  of  Y)  levels  k and  k'  of  Z are 
matched  with  all  pairs  of  adjacent  rows  of  X,  and  log  odds-ratios  of 
$2k'  ” ^2k  w^en  (within  levels  of  X)  levels  k and  k1  are  matched  with 
all  pairs  of  adjacent  columns  of  Y.  This  latter  model  has  linear 
effects  of  the  ordinal  variables  in  all  the  pairwise  associations. 

The  multiplicative  effects  model  (3.27),  however,  can  be  used  with 
nominal  or  ordinal  variables.  For  instance,  if  the  variables  X and  Y 


are  ordinal  the  p-parameters  can  be  interpreted  as  scores,  and  we  will 
be  interested  in  whether  they  are  monotonic.  If  one  set  of  the 
parameters  is  not  ordered,  say  the  scores  for  X,  then  we  can  interpret 
them  as  row-effects  as  if  Xwas  nominal. 

All  the  models  discussed  in  this  section  can  be  further  refined 
to  describe  higher  order  interactions  and  nonuniform  associations. 
However,  in  many  applications  the  proposed  models  will  suffice  for 
describing  the  most  basic  patterns  in  the  data. 


3.7.  Comments  on  Model  Selection 

Several  models,  of  varying  degrees  of  complexity,  have  been  dis- 
cussed in  this  chapter.  Unlike  the  standard  log! inear  models  where  the 
variables  are  treated  as  nominal,  the  user  has  a rich  variety  of  models 
to  choose  from  if  the  ordinal  nature  of  the  variables  is  taken  into 
account.  However,  this  great  variety  of  models  can  make  model  selection 
more  difficult.  In  other  words,  the  researcher  may  be  confused  as  to 
what  strategy  to  take  in  selecting  and  fitting  these  models. 


57 


One  approach  is  to  let  the  goodness-of-fit  standard  (nominal) 
loglinear  models  with  residual  patterns  suggest  ways  of  modelling 
ordinal ity.  For  instance,  suppose  the  standard  loglinear  model 
allowing  all  pairwise  associations  (denoted  by  (X^^X^, . . . ,Xd  lXd)) 
fits  well.  Then,  models  for  ordinal  variables  displaying  no  three 
factor  interaction  can  be  considered.  In  particular,  perhaps  model 
(3.24),  or  model  (3.26)  if  not  all  the  variables  are  ordinal,  or  their 
generalizations  having  homogeneous  additive  or  multiplicative  effects 
will  fit  well.  Such  models  will  be  more  parsimonious  than  the  standard 
models  and  perhaps  easier  to  interpret.  Similarly,  if  a standard  log- 
linear  model  with  certain  interactions  fits  well,  we  may  be  able  to 
obtain  a good-fitting,  simpler  model  by  adding  certain  interaction 
terms  (with  linear  ordinal  effects)  to  the  mentioned  models  above. 

To  illustrate  this  strategy,  we  consider  Table  3.4,  taken  from 
Forthofer  and  Lehnen  (1981,  p.  21).  This  table  concerns  associations 
among  A = age,  S = smoking  status,  and  B = breathing  test  for  Caucasians 
in  certain  industrial  plants  in  Houston,  Texas, during  1974-1975.  The 

standard  loglinear  model  (AS,  AB,  BS)  which  allows  all  bivariate 

2 

associations  fits  poorly,  with  G = 25.9  based  on  df  = 4,  so  we  consider 
models  which  exhibit  three  factor  interaction  terms.  For  integer  scores, 
the  uniform  interaction  model  (3.13)  has  only  one  more  parameter  BABS, 
but  fits  strikingly  better,  as  G2  = 2.7  based  on  df  = 3.  The  ML 
estimate  of  the  uniform  local  interaction  measure  log  0...  is  6ABS  = .83, 

1 J K 

with  standard  error  of  .19.  Hence,  the  association  between  smoking 
status  and  breath  test  is  more  positive  at  the  higher  age  level.  If 
we  are  mainly  interested  in  the  B-S  association  or  the  AB(AS)  marginal 
table  is  fixed  by  sampling  scheme,  we  might  want  to  fit  the 
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heterogeneous,  uniform  B-S  association  model  (i.e.,  integer  scores 

version  of  model  (3.12)).  The  model  does  not  fit  quite  well,  as 
2 

G = 10.8  based  on  df  = 6;  however,  it  yields  the  simple  interpretation 

of  a constant  local  log  odds-ratio  of  1.12  for  the  age  <40  group  and 

2.18  for  the  40-59  age  group.  The  simpler  model  (3.8)  with  integer 

2 

scores  fits  even  worse  since  G =63.1  based  on  df  = 8. 

For  Table  3.1,  the  partial  association  model  (HC,  HP,  CP)  fits 
2 

well,  with  G = 8.1  based  on  df  = 9.  Hence,  we  considered  the  models 

listed  in  Table  3.2,  which  allows  all  bivariate  associations  but 

assumed  no  three-factor  interaction. 

In  many  applications  several  ordinal  models  will  fit  the  data 

reasonably  well.  Of  the  no  three-factor  interaction  models,  it  seems 

to  us  that  the  linear  ordinal -effects  model  (3.27)  is  a natural  first 

choice  for  a researcher.  In  many  (perhaps  most)  applications,  the 

researcher's  theoretical  foundations  will  predict  that  associations 

between  pairs  of  ordinal  variables  are  monotonic,  and  that  levels  of 

nominal  variables  are  stochastically  ordered  with  respect  to  their 

associations  with  ordinal  variables.  Hence,  model  (3.27)  should  often 

provide  a fairly  good  approximation  to  association  patterns,  even  if 

more  complex  models  are  needed  to  pass  formal  goodness-of-fit  tests. 

If  the  model  fits,  particular  associations  can  be  tested  by  the  reduction 
2 

in  G , and  a simpler  model  exhibiting  conditional  independence  between 
certain  pairs  of  variables  may  be  adequate. 

Our  suggestion  of  using  (3.27)  as  a reference  point  for  the 
ordinal  models  is  based  on  considerations  about  the  way  researchers 
commonly  model  quantitative  response  variables  in  regression  models. 
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Unless  research  hypotheses  suggest  otherwise,  the  usual  basic  model  is 
one  in  which  explanatory  variables  have  linear  effects  if  they  are 
quantitative  and  additive  effects  if  they  are  qualitative.  The  re- 
gression coefficients  summarize  the  linear  tendencies  in  associations 
and  interactions  and  are  useful  for  descriptive  purposes  if  we  expect 
the  true  underlying  relationship  to  be  somewhat  more  complex.  It  seems 
to  us  that  models  such  as  (3.27)  and  (3.8)  that  have  linear  effects  in 
the  associations  and  interactions  are  analogous  loglinear  models  to 
the  linear  regression  model  just  described.  Because  of  the  simplicity 
of  interpretation  for  models  with  linear  effects  of  ordinal  variables, 
we  believe  that  it  is  often  nearly  as  informative  to  fit  them  as  to 
fit  more  complex  models  having  additional  effects,  even  if  the  more 
complex  models  are  better  according  to  formal  goodness-of-fit  criteria. 
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Table  3.4.  Cross-classification  of  Houston  industrial  workers  by 
breathing  test  results,  smoking  status,  and  age. 


Breathing  Test  Results 


Age 

Smoking  Status 

Normal 

Borderline 

Abnormal 

<40 

Never  Smoked 

577 

27 

7 

Former  Smoker 

192 

20 

3 

Current  Smoker 

682 

46 

11 

40-59 

Never  Smoked 

164 

4 

0 

Former  Smoker 

145 

15 

7 

Current  Smoker 

245 

47 

27 

Source: 

Forthofer  and  Lehnen  (1981, 

p.  21). 

CHAPTER  FOUR 

TEST  FOR  INDEPENDENCE  IN 
TWO-WAY  TABLES 

To  test  for  independence  of  the  row  and  column  variables  in  an 
rxc  contingency  table,  the  Pearson  chi-square  and  the  likelihood-ratio 
chi-square  statistics  are  often  used.  However,  when  the  row  and 
column  variables  both  have  ordered  categories,  Haberman  (1981)  con- 
sidered tests  based  both  on  the  canonical  correlation  analysis  of 
Fisher  (1940)  and  the  RC  model  (2.5).  He  has  shown  that  the  resulting 
statistics  from  these  two  approaches  are  asymptotically  equivalent 
under  the  hypothesis  of  independence.  Also,  these  test  statistics 
share  the  properties  of  asymptotic  unbiasedness  and  consistency  with 
the  conventional  chi-square  statistics.  Haberman 's  statistics 
have  neither  uniformly  greater  nor  uniformly  smaller  asymptotic  power 
than  the  Pearson  and  likelihood  ratio  chi-square  tests  for  independence. 
He  indicated  that  the  use  of  the  new  statistics  appears  most  appropriate 
when  the  RC  model  holds. 

In  this  chapter,  tests  of  independence  based  on  the  linear-by- 
1 inear  association  model  (2.4)  are  proposed.  The  description  of  these 
statistics  and  their  properties  (section  3)  follows  the  descriptions  of 
the  conventional  chi-square  statistics  and  Haberman 's  statistics 
(sections  1 and  2).  In  section  4,  the  performance  of  the  proposed 
statistics  is  then  compared  with  the  conventional  statistics.  A 
comparison  with  Haberman's  statistics  is  excluded.  This  is  due  to  the 
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fact  that  the  asymptotic  power  of  the  tests  based  on  the  RC  model  can 
be  found  only  in  principle.  Although  an  approximation  is  available  as 
the  number  of  columns  in  a table  goes  to  infinity,  in  this  study  we 
considered  two-way  cross-classifications  only  as  large  as  10x10. 

4.1 . The  Conventional  Statistics 

Consider  an  rxc  contingency  table  with  frequencies  n..,  l<i<r, 

• J 

l^j<c.  Assume  one  of  three  sampling  schemes  (either  Poisson  or 
multinomial  or  product  multinomial).  Since  the  ML  estimate  of  the 
expected  frequencies  are  the  same  under  each  of  the  sampling  schemes 
above,  any  of  the  above  assumptions  about  the  sampling  distribution 
will  lead  to  the  same  results.  Let  N be  the  sample  size  of  the  table. 
Then  the  independence  model  is 


U 


Npi j = exp(y  + X*  + x1.)  for  all  i,j. 


The  usual  tests  of  independence  rely  on  the  Pearson  chi-squares  Xj 


or  the  likelihood-ratio  chi-square  Gj,  where 

*i = z(nij-*i!j)2/Sij 
and 


(4.1a) 


GI  = 2Dlij  lo9  (n^Vnijj).  (4.2a) 

^ij  = ni+n+j/^  are  the  estimated  expected  cell  counts  under  the 
hypothesis  of  independence.  Consider  local  alternatives  H • n.  .-n.  n 
= d-j j/ v^,  where  the  quantities  {d^}  remain  fixed  as  the  sample  size  N 
increases.  It  is  well  known  (Haberman,  1974a,  p.  112)  that  the 
asymptotic  powers  1^  and  of  Xj  and  Gj  respectively,  at  level  a are 
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p(x^(N$?)>x^a).i  = 1.2,  where  v = (r-l)(c-l),  a is  the  upper 

a-point  of  a chi-square  distribution  with  v degrees  of  freedom,  and 
2 2 

Xv(N$i)  is  a noncentral  chi-square  random  variable  with  noncentrality 

2 ? 
parameters  N$..  and  v degrees  of  freedom.  The  {$f}  are  defined  as 


*1  = W+j>  /piAj 


(4.2a) 


and 


$ 


2 

2 


= ^Pi:j (log  Pi:)-log  Pi+-log  p+j)2  . 


(4.2b) 


4.2.  Haberman's  Statistics 

Alternate  tests  of  independence  may  be  based  on  the  RC  model  (2.5) 


mi j = Np.j  = exp  (y  + A*  + Aj  +6  y^j)  . 

Let  {fay}  denote  the  ML  estimates  of  the  expected  frequencies  {m.  .} 

1 =i =r » 1 = J =c , under  the  RC  model.  Then  the  independence  model  may  be 
tested  against  the  RC  model  by  the  Pearson  chi-square  statistic 


2 

XM 


= E( 


m. 


i j 


.-I  \2,~I 

lOjj)  /lUjj 


(4.3a) 


or  the  likelihood-ratio  chi-square  statistic 


(4.3b) 


The  chi-square  test  statistics  just  described  do  not  have  approxi- 
mate chi-square  distributions  under  the  null  hypothesis.  Haberman 

(1981)  derived  the  asymptotic  distribution  of  the  statistics  y2  and  G2 

AM  M 

under  local  alternatives  Ha:  Pij-Pi+P+j  = d.y/N.  Let  H be  the  (rxr) 
matrix  with  elements  h^.,  defined  as  follows: 
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h^i  = l<i,  i'<r  where 

J 

aij  = (pij-pi+p+j)//pi+P+j  ]^c  • 

If  we  let  k = min  (r-l,c-l)  and  p2,  l<j<k  be  the  jth  largest  eigenvalue 
of  H,  then  we  have  $ = trace  (H)  = za?.  = Zp^.  For  local  alternatives, 
Haberman  showed  that  the  asymptotic  p’ower  of  and  G2  at  level  a is 

nM  = P(F(i — 1 ,c-l  ,v%)>x(r-l  ,c-l  ,a)).  (4.4) 


In  the  equation  above,  F(r-1  ,c-l  ,/NgJ  is  the  largest  eigenvalue  of  the 
noncentral  Wishart  matrix  W = {w^  , (r-1  ,c-l ,/%)  lsi<r,  1 < i ' < r } such 


c-1 


that  w. 


n 


Z (v-jj+*/Npj6^)(vi  i j+v/NPj5i  >j),  where  the  {v^}  are 
J 


independent  standard  normal  random  variates,  and  6..  denotes  the 

^ J 

Kronecker  delta.  Critical  values  X(r-l,c-l,a)  of  the  statistic  F des- 
cribed above  can  be  found  in  Table  51  of  Pearson  and  Hartley  (1972). 
Note  that  under  the  hypothesis  of  independence  £ = 0. 

The  joint  distribution  of  the  latent  roots  of  a noncentral  Wishart 


distribution  has  been  derived  by  M.L.  Srivastava  and  E.M.  Carter  (1980). 
However,  only  an  approximation  (pointed  out  by  E.M.  Carter  in  a private 
communication)  of  the  distribution  of  the  largest  root  is  available. 

Let  the  noncentrality  matrix  be  (c-l)A,  and  let  its  roots  be  a^,...^ 
such  that  a-j<a2<. . .<aq  = aq+1  = ...  = ap  = 0.  Let  61  be  the  largest  root 
of  the  noncentral  Wishart  matrix,  and  define  a2  = 2(1+28^.  Then  the 
approximate  distribution  of  ©1  is: 


P[/c-l  (0] -(1+  c_i  P] ) )<z]  - $(z)  + /c-1  [a-j  (l-z2)<j)(z)  + a2<j>(z ) 


+ o(l/c-l )]  as  c + ® 


(4.5) 
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where  a-,  = 4(l+3a,  )/3o3,  a 


and  4>(z),  $(z)  are  the  standard  normal  density  and  distribution  functions, 
respectively. 

The  approximation  is  appropriate  only  when  the  number  of  columns  c 
is  very  large.  We  consider  in  our  power  studies  in  this  chapter  only 
tables  with  at  most  10  columns.  Thus,  we  could  not  make  meaningful 
power  approximations  for  Haberman's  statistics.  (Actually  for  c = 10, 
the  approximation  fails  to  give  accuracy  even  to  the  first  decimal 
place.)  Since  the  fit  of  the  RC  model  requires  many  iterations,  we  did 
not  do  Monte  Carlo  studies  to  obtain  power  approximations  for 
Haberman's  statistics. 

In  summary,  it  is  possible  to  carry  out  tests  of  independence  using 
Haberman's  statistics,  since  their  null  distribution  is  tabulated. 

However,  at  this  time,  power  computations  are  impractical  for  small  c. 

4.3.  Proposed  Statistics 

If  a cross-classification  has  ordered  categories  for  both  variables, 
alternate  tests  of  independence  may  be  based  on  a parsimonious  model 
for  the  table.  For  instance,  simple  tests  can  be  based  on  the  linear- 
by-1 inear  association  model  (2.4),  which  has  one  parameter  more  than 
the  independence  model. 

Given  model  (2.4),  testing  for  independence  is  equivalent  to 
testing  the  following  composite  hypothesis: 

H : 6 = 0 vs.  H : g * 0 

0 d 


parameter.  In  what  follows  through  page  72,  the  parameters,  the  statistics 
and  the  "o"  are  vectors. 


where  6 = (g,A)T,  with  A = (y,A*, 


• • • 9 


9 • • • 9 


A^_i )T  being  a nuisance 
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In  what  follows,  we  will  describe  a general  theory,  based  on  Cox 
and  Hinkley  (1974),  that  will  suggest  three  statistics  for  this  com- 
posite null  hypothesis  that  are  asymptotically  equivalent  under  H . 


In  particular,  we  will  derive  two  statistics  that  are  asymptotically 
equivalent  to  the  likelihood  ratio  statistic  based  on  model  (2.4). 

Suppose  the  random  variable  Y has  p.d.f  f(y;  6, A)  satisfying 
certain  regularity  conditions  (see  Appendix  1).  Suppose  we  wish  to 
test  Hq:  B = Bq  against  Ha:  g * B0>  with  A as  a nuisance  parameter. 
Under  the  null  hypothesis,  A is  estimated  by  maximum  likelihood  with 
B = Bq.  We  use  the  notation  to  denote  the  ML  estimate  of  A subject 
to  e = 30>  and  A to  denote  the  unrestricted  ML  estimate  of  A.  If  we 
have  N i.i.d  random  variables  with  p.d.f  f(y;  g,A),  then  the  loglikeli- 
hood  function  is 


Also  the  efficient  score  for  one  observation,  and  the  total  score  for 
N observations  are  respectively  given  by  the  following  column  vectors: 


In  the  regular  case  under  discussion,  we  have  the  following  properties: 


o 


N 


N 


N 


E(U-j  (e)  ;0)  = E(U.(e);e)  = 0 
E(U1(e)u}(0);e)  = E(-|g  u{(0);0)  = 1^0) 

E(U.  ( 0 ) uT ( 0 ) ;0)  = E (-  |g  uT ( 0 ) ; 0 ) = (0)  = 1.(0) 


(4.7b) 


(4.7a) 


(4.6) 
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The  symbol  1.(0)  denotes  the  information  matrix. 

The  ML  estimate  0 satisfies  the  closely  related  likelihood  equations 

M§)  = 0.  (4.8) 

Now  we  will  derive  the  asymptotic  equivalence  of  three  test 
statistics  for  this  composite  hypothesis  and  local  alternatives.  Let 
W denote  the  likelihood-ratio  statistic  for  testing  H , that  is, 
e5W  = f(6,A;Y)/f(e0,A0;Y).  If  we  expand  W about  the  parameter 
60  = (Bq,A)T,  we  obtain 

W = 2{£(6,A;Y)  -£(B0,Aq;Y)}  = 2{£(b,A;Y)  - £(Bq,A;Y)} 

-2{£(b0,^0;Y)  - £(B0,A;Y)} 


/A  \ 

B-B0 

T 

fl  i"*”  1 

BB  BA 

/A  \ 
B-B0 

’ 0 

T 

fl  IT 

BB  BA 

’ 0 ' 

A 

!ba!aa 

A 

A 

A-A 

A-A 

[Vxj 

.Vaa 

M 

+ eN  * 


(4.9) 


where  the  information  matrix  1.(0  ) = 

0 


(I  iT 

BB  BA 
!ba!aa 


is  partitioned  according 


to  the  partition  (B,A)  of  0.  The  term  e^  includes  powers  and  products 
of  (8-BQ),  (A-A),  and  (Aq-A)  (which  are  Op(l)  by  consistency  of  ML 
estimates)  with  the  third  partial  derivatives  of  the  logl ikel ihood 
function  (which  are  assumed  to  be  bounded  in  the  neighborhood  of  the 
true  value  0Q).  Thus,  e^  is  Op(l).  Equation  (4.9)  reduces  to 

w - (B-B0)TIeB(S-fS0)  + 2(§-60)TI^(X-a)  + (x-a)ti„(a-a) 


AA 


- <H>>  VK>  +V1)- 
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Now  W can  be  simplified  further.  To  do  so,  we  expand  U.(g,A)  and 
u-(B0>y  about  the  Parameter  e0  = (Bq,A)T,  to  find  a relation  between 


(A-A)  and  (Aq-A).  This  gives 


U.(B, A)  = U.(0Q)  + I.(eQ) 


’ A 

B-Bc 

A 

A-A 


+ Opd)  = u.(0o) 


WK>  + 

ISA<B-60)  * 


+ OpC ) . 


and 


MB0.Xq)  = U.(0O)  + I.(0Q) 


f 0 

VA 


+ op(i ) = u.(6o) 


(4.10a) 


XAX(VX) 


+ o (l)  . 
P 


(4.10b) 


Now  the  left  hand  sides  of  equations  (4.10)  are  zero,  by  the  likelihood 
equations  (4.8),  so  that  we  may  equate  one  of  the  components  of  the 
right  hand  side  to  obtain  the  relation: 


!»(*„-*>  - + WS-fU  + °»0) 


AAV  o 
so  that 


BAVM  Ho 1 


r~l 


Xo-X"X-X  + IAiVe-eo)+V1) 


(4.11) 


assuming  IAX  to  be  nonsingular.  Substituting  (4.11)  into  (4.9)  we  obtain 
W = (B-30)  (Ig3"IexIAAI8A^^"eo^  + °p^  * 


(4.12) 
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The  expansion  of  W leads  to  two  asymptotically  equivalent  test  statistics. 
The  first  is  immediate  from  (4.12)  and  is 

We  = (6-60)TI63(6-60),  (4.13) 


where  I33  = 1^  - Hence  1^  is  obtained  from  a partition  of 

1.1 (0O),  the  inverse  of  the  information  matrix  I.(eo).  We  refer  to  the 
statistic  (4.13)  as  the  ML  statistic. 

A second  statistic  equivalent  to  W is  derived  from  the  representa- 
tion of  (B-e0)  in  terms  of  U . ( 6Q » X ) . The  expansion  of  Me  ,A  ) 
about  the  parameter  point  90=(B0>A)T  gives  the  following  equation: 


Mb0,x0>  =U-Ce0)  + fe  uT(s0.x)(e0-e0)  + o 0)  . 


Now  the  left  hand  side  of  the  above  equation  is  zero,  by  the  likelihood 
equations  (4.8),  so  we  may  rewrite  the  above  equation  as 


u-(e0>  = • !e  u-(6o)(Veo)  + °p(1)- 

Multiplying  both  sides  by  l"1(9())  = N_1 1"1  (0O ) we  obtain 

i:1(e0)u.(6o)  . -n"'i^(60)  |e  ul(e0)(§0-e0)  +o  (n'1) 


'[-"'V'Vfe  .f.uI(0o>](V9o>  + “p'N*1)  • 

J * 

Since  E(U1(eo))=0  and  E(u}(e())U1  (0Q)]  = -E[fg  u{(0Q)]  = I]  (eQ),  then 

var[U1(0Q)]  = E[u|(eo)u|(eo)]  - ET[U(0o)]E(U(0o) ) = !■]  (0q)  • From  the 
weak  law  of  large  numbers  applied  to  the  ratio  on  the  right  side  of 

11  r\  N j 

the  equation  above,  we  obtain  that  N“  17  (0  )-2—  z U . (0  ) = 1 + o (1 ) . 

I O 30  j_i  J 0 p'  ' 

Therefore,  the  above  equation  reduces  to: 


70 


(vV'V1 ' I~1  (e0)u- (e0)  • 

In  our  case  we  have  under  the  null  hypothesis  3 = 


f > 

A 

8-8 

f T 

I I 

-1 

' > 
Un 

'jeejSAT 

' " 

u 

0 

— 

88  8A 

8 

+ o (1)  = 

8 

A 

A -A 

T T 

U, 

Pv  1 

jBAjAA 

u 

0 

6A  AA 

A 

UA 

where  U. (0  ) = 
o 


'U 


+ °p{1)> 


is  partitioned  according  to  the  partition  of 


(4.14) 


0 = (6, A)  , and  all  the  quantities  on  the  right  of  (4.14)  are  evaluated 
at  e0  = (B0^)T.  Equation  (4.14)  gives  the  desired  representation  of 

§ - 80  in  terms  of  U.(80,A)  as 


b - e0  - ibbu6  + i6x  ux  + o (i)  , 


(4.15) 


where  the  superscripts  indicate  the  corresponding  submatrices  of  the 

inverse  matrix  IT"*  (6  ) . 

o 

From  (4.12)  and  (4.15)  we  see  that  another  statistic  asymptotically 
equivalent  to  W is  given  by: 


WU(A)  = [I33U6 


(4.16) 


Of  course,  W^(\)  involves  the  unknown  value  A of  the  nuisance  parameter, 
but  the  substitution  of  AQ  for  A in  (4.16)  adds  only  a term  of  order 
OpO)  under  Hq.  This  follows  from  the  fact  that  W (A)  is  a continuous 
function  of  A,  and  that  A , the  ML  estimate  of  A,  is  a consistent 
estimator.  That  is,  AQ  - A = op(l)  implies  W(AQ)  - W(A)  = o (1). 
Therefore,  a usable  form  for  this  statistic  Wu  is  obtained  from  Wu ( A) 
by  replacing  A by  A,  so  that 
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W = W (A  ). 
u uv  o' 


Finally,  since  U.(B  = 0,  W„  simplifies  to 

A 0 O 11 


Wu  * 


(4.17) 


(4.18) 


We  refer  to  this  statistic  as  the  score  statistic. 

By  the  asymptotic  equivalence  of  W,  Wg,  and  W , each  statistic  has 
an  approximate  chi-square  distribution  based  on  one  degree  of  freedom. 
This  follows  from  the  asymptotic  distribution  of  the  likelihood-ratio 
statistic  or  the  asymptotic  normality  of  B and  IL.  The  distribution  is 

p 

central  under  Hq.  For  local  alternatives  BN  = B + 6N,  where  6^  = 5 
and  the  vector  remains  fixed  as  N gets  large,  we  have  for  0^  = (Bn»X)T 


^y(y>6|\j)  fY(y,eo)  ^ 


T sfy0'.9o) 


30. 


+ o(6n) 


(4.19) 


9log  fv(y,0  ) 

Now  E(U.(0O);0N)  =/ — fy(y,eN)  dy, 

and  substituting  (4.19)  in  the  equation  above  we  obtain 
E(u.(e0);eN)  = £(u.(e0),e0)  + «Ji.(e0)  + o(«N). 


The  first  term  on  the  right  hand  side  vanishes  (equation  (4.6)),  so  that 

under  local  alternatives  we  have  E(U.(0j,0j  = fill  - (0  ) + o(5j,  and 

similarly  we  obtain  that  var[U. (0) ;0M]  = E[U.(0  )UT(0  );0j 

O N 0 O N 

-E(U. (0Q) ; 0^] ET[ U . (0Q) ;0N]  = I . ( 0O ) + o ( 6^ ) . Now,  using  the  asymptotic 

normality  of  11.(0),  we  conclude  that,  under  local  alternatives,  the 

score  statistic  W , as  well  as  the  statistics  W and  W , have  an 
^ 6 
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asymptotic  chi-square  distribution  with  noncentrality  parameter 

«NTUe)6N  = 

For  the  linear-by-1 inear  association  model  the  three  statistics 
will  reduce  to  simple  forms.  In  the  following  subsections,  we  will 
derive  the  forms  taken  by  those  statistics  with  their  noncentrality 
parameters  under  local  alternatives. 


4.3.1 . Likelihood  Ratio  Statistic 

In  our  context  the  likelihood  ratio  statistic  is 


W = lostJHj/aJj)  - 2Zn1;j  - ZEn^  logfr,,^) 


» 


where  the  (m^j*  l=i=r,  l^j^c}  are  the  estimated  expected  frequencies 

under  the  1 inear-by-1 inear  association  model,  and  the  {£?.}  are  as 

■ J 

defined  previously.  That  is,  W is  the  likelihood-ratio  chi-square  test 

for  independence  given  that  model  (2.4)  holds.  Assuming  the  model 

holds,  we  know  that  Gy  will  be  asymptotically  distributed  as  a central 

chi-square  with  df  = rc  - r - c.  If  the  alternative  hypothesis  moves 

steadily  closer  to  the  null  hypothesis  (i.e.,  p.  . - p.  p . = d..//N), 

2 J ^ T J 

then  Gj  will  have  an  asymptotic  noncentral  chi-square  distribution  with 

df  = (r-1 ) (c-1 ) and  noncentrality  parameter  from  section  1.  Thus, 
based  on  results  in  Bishop  et  al . (1975,  p.  525)  W will  have  an 
asymptotic  noncentral  chi-square  distribution  with  df  = 1 and  non- 
centrality parameter  N^.  We  note  that  W has  the  same  noncentrality 
2 

as  Gj  does,  but  it  is  focused  here  on  df  = 1 rather  than  df  = (r-1) (c-1). 
Thus,  the  asymptotic  power  of  W at  level  a,  is  given  by: 


(4.20) 


n3  ■ - p(x>22)>4>a) 

2 2 

where  x-|  (N$2)  is  a noncentral  chi-square  random  variable  with  non- 
centrality parameter  N$2  and  df  = 1 , and  x2  a is  the  upper-a  point  of 
a chi-square  distribution  with  df  = 1 . 


4.3.2.  The  ML  Statistic 

For  testing  Hq:  8=0,  the  ML  statistic  simplifies  to 

we  - (8/ag)2, 

where  8 is  the  ML  estimate  of  8>  and  ag  its  estimated  standard  deviation 
under  Hq.  Consider  local  alternatives  8 = 6 //N.  From  model  (2.1)  and 

p 

model  (2.4)  we  note  that 


log  Pij  - log  p.+p+j  = 6(ui-u)(vj-v). 
If  we  set  pij.  - p.+p+j  = d^./ZN,  then 


(4.21) 


log  (Pij/Pi+P+J.)  = log 


u Lj£pi-p-j 


pi+p+o 


= log 


1 + 


d.. 

u 


XNp1+p+j 


d. 


^PT+P+j 


LL-  + o(N'i), 


(4.22) 


From  equations  (4.21)  and  (4.22)  we  remark  that  considering  local 
alternatives  8 = Sg/v'ffi’  is  equivalent  to  considering  local  alternatives 


Pij  " Pi+P+j  = dij/^+  o(N-i),  where  djj  = ^(u.-u) (v ,-v)p.+p  . 


By  the  asymptotic  normality  of  8 (see  Appendix  1),  we  see  that,  under 
local  alternatives,  the  statistic  Wg  will  have  a noncentral  chi-square 
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distribution  with  noncentrality  parameter  N$2,  where  4>2  = (B/o)2  and 
2 _ 2 e e 

o - Nag  . Hence,  the  asymptotic  power  of  Wg,  at  level  a can  be  expressed 

as 

n4  = p(we>xl2,a)  = p(x^(N4>g)>x^a)  • (4.23) 

We  remark  that  n4  can  be  expressed  in  terms  of  the  standard  normal 
distribution 

n4  = p(|6/a'e!>zrj/2)  = P(g/crg>2a/2)  + P(B/ag<-zo/2) 


fS-6fl//N 

> z 

V 

4-  D 

P-v* . **  1 

Og  a/2 

/Nag 

BJ 

• r 

°§  ' ‘a/2 

BJ 

where  $ is  the  standard  normal  distribution  function  and  z its  upper  a 
point.  The  above  representation  in  terms  of  the  standard  normal  dis- 
tribution is  helpful,  since  it  can  be  used  in  an  analogous  way  if  we 
wish  to  perform  only  a one-sided  test. 

4.3.3.  The  Score  Statistic 

The  kernel  of  the  log  likelihood  function  under  the  multinomial 
sampling  scheme  is  log  l_M  = Zn^  log  . subject  to  Emi  . = N.  If  the 
stated  model  is  the  1 inear-by-1 inear  association  model,  the  equation 
above  becomes 

lo9  L|y|  = 2n1.j(y+A^+Aj+S(u-.-u)(Vj-v))  subject  to 
Zexp(y+A^+Aj+3(ui-‘u)(Vj-'v))=  N.  fence. 
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uB(e,x)  = 


Slog  L 


M 
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= J(ui-u)(»rv)(niJ.-m1j). 


Now  under  the  hypothesis  of  independence  (6=0),  we  have 

V°’V  ‘ E<V“)(Vj-7)(nij-ftij)  = I(u.-U)(Vj-v')(njj.-n^+n+j/N) 

= N[r(ui-Tr)(vj-V)(f1j-f1+f+.)] 

' N[ju1vj(fij-f1+f+j)  - “VVV+J>  -v^i(firfi+ftj) 

+ “«<W+J» 


= N[Euivj(f1j-fi+f+j)a, 

where  the  {f^.  = n^./N}  are  the  observed  cell  proportions. 

Let  us  now  consider  how  to  obtain  the  asymptotic  variance  of  this 
type  of  statistic.  Let  T = N_1UQ(0X ) = Eu.v  .(f . .-f . f . = q(f) 

o 6 o'  i j'  ij  i+  yv-' 

say,  where  f_  = (f^  .f^,.  ••  »frc)  is  the  rc-vector  of  observed  cell 
proportions.  Linder  multinomial  sampling  we  have: 

/N  (f-£)  X N(0,z(p)), 

where  £ - (P] i » P] 2* * * ' Ppc^  1 s ^he  *"c“vect°r  of  cell  proportions,  and 
£(£.)  = ^aij}k^,l*i jk£r,l£j,£<c}  is  the  (rcxrc)  covariance  matrix  with 

aij,ij  = pij(1-pij^’  and  aij,k£  = ~pi jpk£  for  i?tk  and/or 

To  derive  the  asymptotic  variance  of  Tq,  we  will  use  the  delta  method. 
That  is, 


*^(g(I)-g(£))  X N(0,dgTz(p)dg), 

8g 

where  dg  - {—  ,l^k£r,lSi!<c}  is  the  rc-column  of  partial  derivatives, 
k £ 


76 


3g_(pJ 

9p^"  = V*  " ukmr  " W 
where  mc  = Ev^p^,  and  mr  = Lu^p^.  Thus 

a2(T0)  = N-1dsTz(2.)dg_  = N-1  [lpij.  (u^-u^-v^)2 

- (EPij(Vj-uiVvjmc))2]*  (4.24) 

Under  the  hypothesis  of  independence  p..  = p.  p .,  (4.24)  reduces  to: 

I vJ  I J 

°o<To>  = N-1[(zufpi+)(i;v?P+j)  - (su1p1+)2(zv|p+j) 

- (lvjP+j)2(Zu2p.+)  + (ZulP.+)2(2VjP<.j)2]  . 

Hence  the  approximate  large  sample  variance  TQ  is 

°o<To>  ■ N'’l>iPi+  - UuiPl+)2][ZV2p+;j  - (EVjP+J)2]. 

involving  the  product  of  the  variances  of  the  marginal  distributions. 
Therefore,  the  score  statistic  can  be  expressed  in  the  form 


W..  = 


a2(T  ) 
o o 


N[EUiVi(f.i-fitf+i)]: 


[KV<Su1fi+>  3 


We  remark  that  one  attractive  feature  of  the  score  statistic  Wu  is  that 
we  can  compute  it  without  fitting  the  model. 

Now  we  will  obtain  the  noncentrality  parameter  for  the  score 
statistic.  In  its  normal  form  it  can  be  expressed  as: 


To 

W 


ayj^ij-pjj1  _ ya  °,<v  , 2v.j(Pi.i-p°i> 


a TTj 

o o' 


°o<U 


w 
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where  the  {p0..}  and  {pa..}  are  the  true  cell  proportions  under  the  null 
hypothesis  and  alternative  hypothesis  respectively,  T = Zu.v.(f . .-p® .)» 
and  ct  (T  ) is  a function  of  the  {p?.}  as  in  (4.24).  Consider  local 

da  I J 

alternatives  p®.  - p?.  = d ../At.  The  quantity  a (T  )/a  (T  ) tends  to  1 

i J i J I J a a 0 0 

as  N becomes  large,  and  Ta/oa(Ta)  tends  to  become  normally  distributed 

with  zero  mean  and  unit  variance.  Consequently,  so  does  T /a0(T  ) so 

that  T /c  (T  ) becomes  normally  distributed  with  mean  Zu.v.(p® .-p?.)/a  (T  ) 
ooo  k-ij"  ox  o' 

and  unit  variance.  Thus,  under  local  alternatives  [T  /a  (T  )12  is 

0 0 0 

noncentral  chi-square  distributed  with  df  = 1 and  noncentrality 
2 

parameter  N$u  where 


However,  for  the  score  statistic  Wu  the  cell  proportions  are 
estimated  from  the  sample.  Let  L(p)  be  the  likelihood  ratio  statistic 
for  testing  independence  when  the  cell  proportions  are  assumed  to  be 
known,  and  L(p)  the  likelihood  ratio  statistic  for  testing  independence 
when  the  cell  proportions  are  estimated  from  the  sample.  Wald  (1943) 
has  shown,  under  local  alternatives,  both  statistics  will  have  the  same 
distribution  as  the  sample  size  goes  to  infinity.  Using  the  asymptotic 
equivalence  of  W and  Wg,  we  conclude  that  Wu  is  noncentral  chi-square 

O 

with  df  = 1 , and  noncentrality  parameter  Therefore,  the  asymptotic 

power  of  Wu,  at  level  a is 

J5  ■ P<Wu<a>  = P<x(<N*u>>Xl,a>- 


(4.25) 
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Since  the  tail  probability  of  a noncentral  chi-square  random 

variable  is  an  increasing  function  of  the  noncentrality  parameter 

(Johnson  and  Kotz,  1970,  p.  135),  all  tests  under  study  are  asymptotically 

unbiased.  Consider  fixed  alternatives  H, : 8=6  * 0.  From  the 

a a 

asymptotic  normality  of  § the  asymptotic  variance  of  6 has  the  form 
2 2 

ag  = c /N.  For  the  ML  statistic  Wg,  we  have  for  large  N 


p(|g/aa|>2a /Z)  = PCt/ofZa/z) 


P(B/ag<- 


Za/2  } 


= P 


M<z  - £- 

ag  a/2  ag 


+ P 


ag  a/2  ag 


1 -$ 


z - £- 

a/ 2 a 


B 


+ $ 


-Z  - £-1 
a/ 2 ag 


+ o(n~*) 


= 1 -$ 


za/2"  -^1  + * 


o//N 


0/2  o//N 


+ o(n~h' 


Hence,  as  N -*■  ®,  P(  | B/ag|>za/,J-*-l,  so  that  the  ML  statistic  is  asymptotically 
consistent.  It  can  be  also  shown  that  both  statistics  W and  Wu  are 
asymptotically  consistent.  Therefore,  all  test  statistics  under  investi- 
gation in  this  chapter  are  asymptotical ly  unbiased  and  asymptotically  con- 
sistent. 


4.4.  Example 

The  4x3  cross-classification  (Grizzle  et  al . , 1969)  in  Table  4.1 
gives  the  dumping  severity  for  four  ulcer  operations.  The  values  in 
parentheses  are  the  ML  estimates  of  expected  counts  under  the  uniform 
association  model.  As  mentioned  earlier,  we  do  not  have  to  fit  the 
model  to  compute  the  score  statistic.  However,  the  fit  of  model  (2.4) 
is  necessary  to  compute  the  ML  statistic,  and  in  addition  to  that  the 
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independence  model  is  fitted  to  compute  the  likelihood-ratio  statistic. 

The  usual  likelihood  chi-square  statistic  for  independence  has  its 
2 

value  Gj  = 10.88  based  on  df  = 6,  and  has  an  observed  significance 

level  .092.  The  statistic  G2  fails  to  reject  the  hypothesis  of 

independence  at  a = .05.  The  RC  model  fits  the  data  quite  well  with 

Grc  = 2.85  based  on  df  = 2.  Haberman's  statistic  G( I | RC)  = G?  - gL  = 8.03, 

1 KL 

and  its  critical  values  X(r-l,c-l,a)  found  from  Table  51  (Pearson  and 
Hartley  (1972))  are  X(3,2;05)  = 10.74  and  X(3,2,.01)  = 14.57.  Thus, 
Haberman's  statistic  fails  to  reject  the  hypothesis  of  independence  at 
a = .05.  This  agrees  with  the  conclusion  based  on  G2. 

Consider  now  the  proposed  statistics.  The  uniform  association  model 
provides  a good  fit  with  G^  = 4.59  based  on  df  = 5.  The  modest  improve- 
ment in  model  fit  with  the  RC  model  does  not  outweigh  its  additional 
complexity  over  the  uniform  association  model.  The  estimated  constant 
value  of  the  local  log  odds  is  found  to  equal  .163,  and  its  estimated 
standard  deviation  equals  .065.  The  values  of  the  proposed  statistics 

are  w = G2  - G2  = 6.29,  Wg  = (g/Sg)2  = 6.29,  and  Wu  = (TQ/aT  )2  = 6.23. 

o 

The  chi-square  statistics  W,  Wg,  and  Wu  based  on  a single  degree  of 
freedom  are  approximately  equal  with  observed  significance  level  .01. 

They  suggest  strongly  that  we  can  reject  the  hypothesis  of  independence. 

So,  for  the  sample  table,  the  conventional  chi-square  statistic  G2  for 
testing  independence  and  Haberman's  statistic  G2 ( I | RC ) , treating  the 
categorical  variables  as  if  they  were  nominal  in  scale,  do  not  seem  to 
present  strong  evidence  against  the  hypothesis  of  independence.  But  the 
proposed  statistics,  using  the  natural  order  of  the  categories,  suggest 
strongly  that  there  is  a trend  in  the  data. 
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4.5.  Power  Comparison 

In  this  section  we  conduct  a numerical  power  study  to  investigate 
the  performance  of  our  proposed  statistics,  comparing  them  to  the  con- 
ventional chi-square  statistics.  Under  the  hypothesis  of  independence 
and  for  local  alternatives,  the  conventional  statistics  are  asymptotically 
equivalent.  The  proposed  statistics  are  also  asymptotically  equivalent 
under  these  hypotheses,  but  not  asymptotically  equivalent  to  the  con- 
ventional statistics.  Since  the  former  statistics  appear  to  perform 
well  when  the  categorical  variables  are  nearly  independent  (Haberman, 
1981),  we  consider  two-way  tables  where  the  variables  manifest  a weak 
association.  The  study  of  the  power  comparison  is  numerical  through  the 
use  of  noncentrality  formulae  of  Section  3.  These  are  approximate 
powers,  since  the  noncentrality  formulae  are  valid  only  under  local 
alternatives  and  for  large  sample  sizes. 

4.5.1 . Description 

Two-way  tables  with  three  different  sizes  (3x3,  5x5,  and  10x10) 
are  investigated.  The  investigation  includes  the  computation  of  the 
different  asymptotic  powers  {ni : 1 <i<5}  for  three  different  levels 
(.01,  .05,  and  .10),  and  for  different  sample  sizes  ranging  from  50  to 
50,000. 

The  formula  for  the  standard  bivariate  normal  density  function 
(i.e.,  y^  = piy  = 0,  and  = Oy  = 1 ) at  X = x and  Y = y includes  the 
product  xyp/(l-p  ),  where  p is  the  coefficient  of  correlation.  The 
corresponding  products  in  the  exponent  of  the  1 inear-by-1 inear  associa- 
tion model  and  the  RC  model  are  (3u.v.  and  gy.v.  respectively. 

Thus,  due  to  this  similarity,  we  consider  an  underlying  standard  bivariate 


normal  distribution  with  a small  correlation  coefficient  (.05,. 20) 
so  that  the  variables  are  weakly  associated. 
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Treating  both  classifications  the  same  way,  we  discretize  the 
bivariate  distribution  according  to  two  different  spacings  of  the 
category  boundaries.  First,  we  consider  an  "equal"  spacing  of  the 
category  boundaries.  We  refer  to  the  "equal"  spacing  as  case  A. 

For  instance,  for  a 5x5  table  we  mean  by  equal  spacing  the  discretization 
of  the  bivariate  distribution  at  the  points  X = -1 .5,-. 5, .5,1 .5  and 
Y = -1.5,  -.5, .5, 1.5.  Then  the  cell  proportions  of  the  table  under 
consideration  are  as  follows:  P11=P(X<-1 . 5, Y<-1 .5),  P12=P(X<-1 .5,-1 .5<Y<-0.5), 

P13=P(X<-1.5,-0.5<Y<0.5),  p14=P(X<-1.5,0.5<Y<1.5),  P15=P(X<-1 .5,Y>1 .5), 
P21=P(-1.5<X<-.5,Y<-1.5),  p22=P ( -1 . 5<X<- .5,-1 .5<Y<-.5), 
P23=P(-1.5<X<-.5,-.5<Y<.5),  p24=P ( -1 . 5<X<-. 5, . 5<Y<1 . 5) , 

P25=p(_1 •5<X<-.5,Y>1 .5),  p31=P(-.5<X<.5,Y<-1.5),  P32=p ( " - 5<X< . 5 , -1 .5<Y<-.5), 
P33=P("‘ 5<X<* 5,“* 5<Y<* 5) , P34=p(--5<X<.5,.5<Y<1.5), 

p35=p(-* 5<x<- 5>Y>1 *5)»  P41=P(.5<X<1.5,Y<-1.5),  p42=P(.5<X<1 .5,-1 .5<Y<-.5), 
P43=P(.5<X<1.5,-.5<Y<.5),  P44=P ( . 5<X<1 .5, . 5<Y<1 .5), 

P45=p(.5<X<1.5,Y>1.5),  p51=P(X>1.5,Y<-1.5),  P52=P(X>1 . 5,-1 . 5<Y<-. 5) , 

P53=P(X>1 •5»_'5<Y<*5)>  P54=p(X>1.5,  .5<Y<1 .5),  p55=P(X>l .5,Y>1 .5) . 

Second,  we  consider  case  B,  where  the  category  boundaries  are  "unequally" 
spaced.  The  cut-points  of  the  category  boundaries  of  the  bivariate 
distribution  for  the  two  spacings  are  displayed  in  Table  4.2.  For  each 
different  spacing,  we  consider  two  underlying  bivariate  distributions. 

That  is,  the  coefficient  of  correlation  p between  the  two  variables 
equals  .05  for  one  case,  and  equals  .20  for  the  other  case.  To  deter- 
mine the  adequacy  of  the  uniform  association  model,  we  fit  model  (2.4) 
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to  the  population  tables  under  investigation.  The  cell  proportions  in 
the  population  tables  are  multiplied  by  1000,  and  the  statistic  Gy  is 
used  to  measure  how  closely  they  agree.  The  results  are  reported  in 
Table  4.3,  and  we  see  that  the  spacing  in  case  A provides  an  excellent 
fit  of  model  (2.4),  whereas  the  spacing  in  case  B gives  a poor  fit. 

Finally,  we  consider  case  C,  in  which  the  cell  proportions  used  for 
the  analysis  are  smoothed  versions  that  satisfy  either  the  uniform 
association  model  or  the  RC  model.  In  the  first  set  of  tables,  the 
uniform  association  model  fits  perfectly.  In  the  second  set  the  RC 
model  fits  exactly,  but  the  uniform  association  model  fits  poorly. 

The  latter  set  enables  us  to  take  into  account  the  case  when  the 
categories  are  not  equally  spaced,  but  the  user  instead  assigns  integer 
scores  to  the  levels  of  the  categories  to  fit  the  model.  In  this  case 
the  RC  model  holds  perfectly  when  for  both  sets  of  categories  the  scores 
(-2,-1 . 5,1  ),(-3,-l  ,0, . 5, . 3) , and  (-3,-2. 5,-1  ,-.5,0, .3,  .8,1 .2,2. 5,3)  are 
assigned  to  the  levels  of  the  3x3,  5x5,  and  10x10  cross-classifications, 
respectively.  For  both  models  the  association  parameter  g is  comparable 
in  value  to  a correlation  coefficient  p of  order  .01  in  the  bivariate 
normal  distribution.  Due  to  the  similarity,  mentioned  earlier,  of  the 
models  and  the  bivariate  normal  distribution,  we  choose  B so  it 
corresponds  to  p = .01  when  equated  in  B = p/(l-p^). 

4.5.2.  Results 

4. 5. 2.1.  Case  A 

The  powers  are  computed  by  using  the  noncentrality  formulae  of  the 
previous  section  and  hence  are  approximate,  moreso  as  the  coefficient 
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of  correlation  p departs  further  from  0.  All  the  computed  powers  in- 
crease as  the  sample  size  N and/or  the  association  p between  the 
categorical  variables  increase.  However,  the  computed  powers  n-j  and 
II2  of  the  conventional  statistics  decrease  as  the  size  of  the  table  gets 
larger,  whereas  the  computed  powers  n3,  n4>  and  n5  for  the  proposed 

statistics  increase.  This  follows  from  the  fact  that  the  noncentrality 
2 

parameters  {4^;  i = 1,2,..., 5}  do  not  vary  much  from  one  table  to  the 
corresponding  larger  table,  and  hence  increasing  degrees  of  freedom 
produce  lower  powers.  That  is,  as  the  size  of  the  table  gets  larger, 
the  proposed  statistics  (based  on  only  df=l ) should  perform  better  than 
the  conventional  statistics,  for  which  df  = (r-l)(c-l)  increases. 

Since  the  noncentrality  parameters  are  practically  equal  in  each  table  of 
fixed  size,  the  proposed  statistics  perform  equally.  The  statistics  W,  We, 
and  Wu  dominate  uniformly  the  statistics  Xj  and  Gj  in  this  present  case. 

For  small  N the  computed  powers  for  the  proposed  statistics  are  minimally 
better  for  the  case  in  which  the  variables  are  weakly  associated 
(p  = -05),  but  when  p = .20,  these  powers  become  substantially  better. 

Since  all  statistics  under  investigation  are  asymptotically 
consistent,  their  asymptotic  powers  tend  to  1 as  we  depart  from  in- 
dependence. To  examine  the  behavior  of  the  proposed  statistics  as  a 
function  of  p,  we  compute  the  relative  type  II  error  rate  Ejj  for  the 
different  values  of  p.  Since  n-|,  and  ji2  are  practically  equal,  and  so 
are  n^,  n^,  and  JI3,we  chose  = (l-n2)/(l-n4).  Other  ratios  can  be 
chosen  as  well.  We  see  from  Table  4.4  that  Ejj  increases  with  p.  This 
suggests  that  the  proposed  statistics  are  more  powerful  for  detecting 
association  as  we  depart  from  independence.  For  all  statistics  in 
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this  study,  the  asymptotic  powers  tend  to  a as  p tends  to  0.  There- 
fore, the  relative  type  II  error  En  tends  to  1 as  p goes  to  0.  For 
specifics  on  power  comparison  we  refer  the  reader  to  Tables  4.5a-4.5f. 

In  summary,  if  there  is  a trend  in  the  data,  the  proposed  statistics  are 
more  powerful  for  detecting  it. 

Since  the  computed  powers  n-j  and  1^  are  relatively  the  same  for 

the  remaining  cases  and  likewise  the  computed  powers  nQ,  II,,  and  nc 

b 

do  not  differ  appreciably  in  most  of  those  cases,  we  will,  for  economy 
sake,  choose  one  statistic  from  each  set,  say  Gj  and  Wg,  to  compare 
powers. 

4. 5.2.2.  Case  B 

For  the  case  that  the  correlation  between  the  categorical  variables 
is  equal  to  .05,  the  conventional  statistics  dominate  the  proposed 
statistics.  This  happens  because  the  uniform  association  model  fits 
the  table  poorly.  The  statistic  Wg  appears  to  perform  better  for 
larger  tables  (10x10)  than  for  smaller  tables  (3x3  and  5x5),  whereas 
the  statistic  Gj  performs  best  when  the  size  of  the  table  is  moderate 
(5x5).  This  might  indicate  that  the  spacing  is  more  severe  for  that 
table  (i.e.,  the  trend  is  less  apparent  in  that  table).  However,  the 
fact  that  the  performance  of  the  statistic  Gj  is  better  for  moderate  size 
tables  than  for  small  tables  may  not  hold  in  general. 

We  note  that  from  Table  4.3  the  model  fits  better  here  than  when 
p = .05,  so  that  the  proposed  statistics  dominate.  This  seems  to 
signal  that  statistics  based  on  model  (2.4)  perform  well  whenever 
there  is  a trend  in  the  table  even  though  the  categories  are  not 
equally  spaced.  It  is  known  that  the  conventional  statistics  per- 
form well  as  we  are  close  to  independence.  If  there  is  a substantial 
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departure  from  independence,  the  proposed  statistics  perform  better 
than  the  conventional  statistics  provided  there  is  a certain  trend  in 
the  data.  Again  we  see  from  Table  4.5b  that  the  performance  of  the 
statistic  Wg  improves  with  the  size  of  the  table,  whereas  the  opposite 
holds  for  the  statistic  Gj.  We  finally  note  that  the  conventional 
statistics  perform  better  in  this  case  than  the  previous  case.  The 
uniform  association  model  is  less  appropriate  in  this  case  than  case  A, 
and  hence  the  proposed  statistics  performed  better  in  the  previous  case. 

4. 5. 2. 3.  Case  C 

First,  we  consider  the  case  in  which  the  uniform  association  model 

holds.  Since  the  association  parameter  B = .01  is  close  to  H , the 

o 

proposed  statistics  are  asymptotically  equivalent,  and  so  are  the  con- 
ventional statistics.  Also  the  statistic  W has  the  same  noncentrality 
2 

parameter  as  Gj  does.  Hence,  the  three  statistics  W,  W , and  W should 

2 2 

perform  better  than  Gj  and  Xj  because  they  are  based  on  fewer  degrees 
of  freedom.  We  see  from  Table  4.7  that  the  computed  power  n4  is 
insignificantly  larger  than  n2  for  the  smallest  sample  size.  The  difference 
between  n2  and  n4  becomes  substantial  as  the  size  of  the  sample  increases, 
particularly  for  large  tables. 

The  case  where  the  RC  model  holds  perfectly  is  considered  to 
complement  case  B.  That  is,  for  moderate  association  in  case  B the 
statistic  We  did  perform  better  than  Gj  even  though  the  categories  are 
somewhat  severely  unequally  spaced.  In  this  present  case  we  consider 
weak  association  ( p=. 01 ) with  milder  unequal  spacing  between  categories. 
Again  the  statistic  Wg  performs  better  than  Gj.  We  see  from  Table  4.8 

r\ 

that  for  both  statistics  Wg  and  Gj,  the  performance  is  better  for 
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moderate  size  of  the  table.  Again  the  latter  result  does  not  hold  in 
general . 

4.5.3.  Conclusions 

For  almost  all  cases  considered  in  this  study,  the  statistics  W, 

Wg,  and  Wu  performed  better  than  the  more  familiar  statistics  Xj  and 
2 

Gj  for  testing  independence  for  two-way  cross-classifications.  The 
performance  improves  as  the  size  of  the  table  gets  larger  and  the 
association  between  the  categorical  'variables  increases. 

In  general,  we  would  expect  the  proposed  statistics  to  perform 
well  whenever  there  is  an  apparent  trend  in  the  data,  provided  the 
category  boundaries  are  not  severely  unequally  spaced.  It  appears 
that  the  conventional  statistics  are  not  as  powerful  for  detecting 
slight  trends  in  the  data,  particularly  for  tables  of  large  size. 
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Table  4.1.  Cross-classification  of  duodenal  ulcer  patients  according 
to  operation  and  dumping  severity. 


Operation 

None 

Dumping  severity 
SI ight 

Moderate 

A 

61 

28 

7 

(62.5) 

(26.2) 

(7.3) 

B 

68 

23 

13 

(62.9) 

(30.9) 

(10.2) 

C 

58 

40 

12 

(61.0) 

(35.3) 

(13.7) 

D 

53 

38 

16 

(53.7) 

(36.6) 

(16.7) 

The  parenthesized  values  are  estimated  expected  counts  under  the 
uniform  association  model. 


Table  4.2.  Location  of  category  boundaries  for  forming  tables  from 
bivariate  normal  distribution. 


3x3 

Equal  spacing 
table  sizes 

5x5 

10x10 

Unequal  spacing 
table  sizes 

3x3  5x5 

10x10 

-2.0 

-2.1 

-1.5 

-1.5 

-1.8 

-1.8 

-1.0 

-1.5 

-1.5 

-0.5 

-0.5 

-0.5 

-1.2 

-1.2 

0.0 

-0.8 

0.5 

0.5 

0.5 

-0.3 

-0.3 

1.0 

0.3 

1.5 

1.5 

.5  1.0 

1.0 

2.0 

1.8 
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Table  4.3.  The  distance  between  the  uniform  association  model  and  the 
cell  proportions  considered  in  this  study. 


Sample  size 
N = 1000 

3x3 

Table  sizes 
5x5 

10x10 

Case  A 

P 

= .05 

.000 

.000 

.000 

P 

= .20 

.000 

.000 

.000 

Case  B 

P 

= .05 

95.900 

383.800 

1164.200 

P 

= .20 

55.000 

215.600 

639.700 

Case  C 

0 

0 

0 

Table  4 

.4. 

Computed  relative  type  II 

error  rates 

Ejj  for  a = .05. 

Table  sizes 

3x3 

5x5 

10x10 

Sample 

size 

p = .05 

p = .20  p = .05 

p = .20 

p = .05  p = .20 

50 

1.006 

1.105  1.011 

1.220 

1.029  1.308 

100 

1.012 

1.237  1.021 

1.559 

1.037  1.847 

200 

1.023 

1.578  1.044 

2.743 

1.057  4.240 

500 

1.062 

3.191  1.121 

17.762 

1.161  86.389 

1000 

1.134 

9.849  1.283 

357.863 

1.392  20818.000 
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Table  4. 5. a.  Computed  powers  when  the  category  boundaries  of  a 3x3' 
table  are  equally  spaced  (p  = .05). 


Noncentrality 

parameters 

*?=.0016 

$2=.0016 

$3=.0016 

$4=. 001 6 

$3=. 0016 

a N 

nl 

n2 

n3 

n4 

n5 

50 

.011 

.011 

.013 

.013 

.013 

100 

.012 

.012 

.016 

.016 

.016 

.01  200 

.015 

.015 

.023 

.023 

.023 

500 

.024 

.024 

.047 

.047 

.047 

1000 

.045 

.045 

.097 

.096 

.096 

50 

.054 

.054 

.059 

.059 

.059 

100 

.058 

.058 

.069 

.069 

.069 

.05  200 

.066 

.066 

.088 

.088 

.088 

500 

.094 

.094 

.147 

.147 

.146 

1000 

.145 

.145 

.246 

.246 

.244 

50 

.106 

.106 

.114 

.114 

.114 

100 

.113 

.113 

.127 

.127 

.127 

.10  200 

.125 

.125 

.155 

.154 

.154 

500 

.165 

.165 

.234 

.233 

.233 

1000 

.235 

.235 

.357 

.356 

.355 

Note: 

Computed  Powers 

nl 

n2 

n3 

n4 

% 

Corresponding 

2 

P 2 

W 

We 

statistics 

XI 

6I 

Wu 
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Table  4.5.b.  Computed  powers  when  the  category  boundaries  of  a 5x5 
table  are  equally  spaced  (p  = .05). 


Noncentrality 

parameters 

$*=.0021 

$2=. 0021 

$3=.0021 

CsJ 

O 

O 

• 

II 

OJ  <3- 

e* 

$g=. 0020 

a N 

ni 

n2 

"3 

n4 

n5 

50 

.010 

.010 

.014 

.014 

.014 

100 

.011 

.011 

.018 

.018 

.018 

.01  200 

.013 

.013 

.027 

.027 

.027 

500 

.017 

.017 

.060 

.060 

.060 

1000 

.028 

.028 

.128 

.128 

.128 

50 

.052 

.052 

.062 

.062 

.062 

100 

.054 

.054 

.074 

.074 

.074 

.05  200 

.059 

.059 

.099 

.099 

.099 

500 

.075 

.075 

.175 

.175 

.175 

1000 

.105 

.105 

.302 

.303 

.300 

50 

.104 

.104 

.117 

.117 

.117 

100 

.107 

.107 

.135 

.135 

.135 

.10  200 

.115 

.115 

.170 

.170 

.169 

500 

.139 

.139 

.269 

.270 

.268 

1000 

.184 

.184 

.420 

.421 

.420 
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Table  4.5.c.  Computed  powers  when  the  category  boundaries  of  a 10x10 
table  are  equally  spaced  (p  = .05). 


Noncentral ity 
parameters 

$*=.0028 

$2=. 0027 

4>3=.0027 

$4=. 0026 

$3=. 0026 

a N 

ni 

n2 

"3 

n4 

n5 

50 

.010 

.010 

.014 

.014 

.014 

100 

.011 

.011 

.019 

.019 

.019 

.01  200 

.011 

.011 

.030 

.030 

.030 

500 

.013 

.013 

.068 

.068 

.066 

1000 

.017 

.017 

.149 

.148 

.148 

50 

.051 

.051 

.064 

.063 

.063 

100 

.052 

.052 

.078 

.077 

.077 

.05  200 

.054 

.054 

.106 

.105 

.105 

500 

.061 

.061 

.192 

.191 

.191 

1000 

.074 

.074 

.336 

.334 

.331 

50 

.102 

.102 

.120 

.120 

.120 

100 

.103 

.103 

.146 

.139 

.139 

.10  200 

.107 

.107 

.179 

.178 

.178 

500 

.118 

.118 

.291 

.290 

.290 

1000 

.138 

.138 

.458 

.455 

.451 
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Table  4.5.d. 

Computed  powers  when  the  category  boundaries 
table  are  equally  spaced  (p  = .20). 

of  a 3x3 

Noncentral ity 
parameters 

**=.0261 

$2=.0260 

$g=.0260 

$4=. 0260 

$5=.0260 

a N 

nl 

n2 

n3 

n4 

n5 

50 

.036 

.036 

.075 

.076 

.076 

100 

.078 

.078 

.167 

.169 

.169 

.01  200 

.200 

.202 

.383 

.386 

.390 

500 

.648 

.651 

.847 

.851 

.851 

1000 

.964 

.965 

.994 

.994 

.994 

50 

.124 

.125 

.207 

.208 

.209 

100 

.214 

.215 

.364 

.366 

.369 

.05  200 

.411 

.413 

.624 

.628 

.630 

500 

.837 

.839 

.950 

.951 

.952 

1000 

.992 

.992 

.9991 

.9992 

.9996 

50 

.208 

.208 

.309 

.310 

.312 

100 

.323 

.324 

.487 

.489 

.490 

.10  200 

.539 

.541 

.736 

.740 

.740 

500 

.902 

.904 

.975 

.976 

.977 

1000 

.997 

.997 

.9997 

.9997 

.9999 
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Table  4.5.e.  Computed  powers  when  the  category  boundaries  of  a 5x5 
table  are  equally  spaced  (p  = .20). 


Noncentrality 

parameters 

$^=.0334 

$2=.0342 

$3=. 0342 

$4=. 0352 

$3=. 0379 

a N 

ni 

n2 

H3 

n4 

n5 

50 

.023 

.024 

.102 

.111 

.109 

100 

.045 

.046 

.233 

.255 

.245 

.01  200 

.117 

.121 

.515 

.553 

.535 

500 

.493 

.506 

.940 

.956 

• .949 

1000 

.925 

.932 

.9995 

.9998 

.9997 

50 

.093 

.094 

.257 

.273 

.266 

100 

.150 

.151 

.456 

.484 

.468 

.05  200 

.290 

.296 

.743 

.774 

.753 

500 

.724 

.735 

.985 

.990 

.989 

1000 

.979 

.982 

.9999 

.998 

.9999 

50 

.166 

.167 

.369 

.387 

.372 

100 

.243 

.247 

.581 

.607 

.590 

.10  200 

.414 

.421 

.834 

.857 

.850 

500 

.821 

.830 

.993 

.9998 

.999 

1000 

.992 

.992 

.9999 

1.000 

.999 
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Table  4.5.f.  Computed  powers  when  the  category  boundaries  of  a 10x10 
table  are  equally  spaced  (p  = .20). 


Noncentrality 

parameters 

^=.0427 

$2=. 0448 

$3=.0448 

$4=. 0452 

CO 

00 

o 

II 

CM  LO 

•& 

a N 

ni 

n2 

n3 

n4 

n5 

50 

.015 

.015 

.111 

.120 

.122 

100 

.023 

.023 

.254 

.276 

.280 

.01  200 

.047 

.046 

.551 

.590 

.591 

500 

.189 

.198 

.956 

.968 

.979 

1000 

.612 

.633 

.9997 

.9998 

.9999 

50 

.069 

.069 

.272 

.288 

.293 

100 

.092 

.093 

.481 

.509 

.519 

.05  200 

.151 

.154 

.772 

.800 

.810 

500 

.408 

.420 

.990 

.993 

.9998 

1000 

.819 

.833 

.9996 

.9999 

1.000 

50 

.130 

.131 

.386 

.405 

.428 

100 

.166 

.168 

.606 

.632 

.639 

.10  200 

.248 

.253 

.856 

.877 

.880 

500 

.544 

.557 

.996 

.997 

.9999 

1000 

.894 

.904 

.9999 

.9999 

1.000 
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Table  4.6. a. 

Computed 
spaced  (p 

powers 
= .05) 

when  the  category  boundaries  are 
• 

unequal ly 

Table  size 

3x3 

5x5 

10x10 

a N 

n2 

n4 

n2 

"4 

n2 

n4 

100 

.012 

.011 

.014 

.011 

.013 

.012 

1000 

.037 

.027 

.118 

.026 

.069 

.061 

.01  2000 

.081 

.048 

.359 

.045 

.224 

.123 

5000 

.286 

.128 

.928 

.118 

.852 

.368 

10000 

.665 

.295 

.9998 

.272 

.9996 

.728 

100 

.056 

.054 

.066 

.054 

.060 

.061 

1000 

.126 

.099 

.292 

.095 

.205 

.168 

.05  2000 

.219 

.149 

.602 

.142 

.455 

.294 

5000 

.517 

.302 

.980 

.286 

.952 

.368 

10000 

.848 

.531 

.9999 

.504 

.9999 

.886 

100 

.110 

.106 

.125 

.106 

.116 

.117 

1000 

.211 

.169 

.417 

.165 

.318 

.261 

.10  2000 

.329 

.237 

.721 

.228 

.592 

.412 

5000 

.643 

.420 

.991 

.402 

.977 

.724 

10000 

.910 

.653 

.9999 

.627 

.9999 

.936 
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Table  4.6.b. 

Computed  powers  when  the  category  boundaries  are  unequally 
spaced  (p  = .20). 

Table  size 

3x3 

5x5 

10x10 

a N 

n2 

n4 

n2 

n4 

n2 

n4 

100 

.063 

.130 

.045 

.205 

.025 

.236 

200 

.157 

.300 

.116 

.461 

.052 

.521 

.01  500 

.539 

.748 

.488 

.910 

.233 

.943 

1000 

.919 

.978 

.923 

.998 

.704 

.9995 

2000 

.9995 

.9999 

.9998 

.9999 

.995 

1.000 

100 

.184 

.306 

.148 

.418 

.098 

.460 

200 

.348 

.537 

.288 

.697 

.167 

.748 

.05  500 

.759 

.900 

.720 

.974 

.466 

.986 

1000 

.977 

.996 

.978 

.9998 

.877 

.9999 

2000 

.9999 

.9999 

.9999 

1.000 

.9993 

1.000 

100 

.286 

.424 

.242 

.547 

.175 

.585 

200 

.474 

.658 

.412 

.797 

.474 

.658 

.10  500 

.846 

.945 

.818 

.980 

.994 

.602 

1000 

.989 

.998 

.990 

.9999 

.933 

.9999 

2000 

.9999 

.9999 

.9999 

1.000 

.9997 

1.000 

97 


Table  4.7. 

Computed 

powers  when 

the  uniform  association  model 

holds. 

Table  size 

3x3 

5x5 

10x10 

a N 

n2 

"4 

"2 

n4 

n2 

n4 

1000 

.012 

.014 

.011 

.014 

.010 

.014 

5000 

.018 

.031 

.013 

.031 

.011 

.032 

.01  10000 

.029 

.057 

.017 

.058 

.013 

.060 

20000 

.057 

.123 

.027 

.124 

.016 

.128 

50000 

.190 

.367 

.076 

.369 

.029 

.382 

1000 

.055 

.062 

.052 

.062 

.051 

.062 

5000 

.076 

.109 

.061 

.109 

.055 

.111 

.05  10000 

.105 

.169 

.074 

.170 

.060 

.174 

20000 

.171 

.293 

.103 

.295 

.071 

.302 

50000 

.396 

.607 

.215 

.612 

.111 

.624 

1000 

.108 

.117 

.103 

.117 

.102 

.118 

5000 

.140 

.182 

.118 

.183 

.108 

.186 

.10  10000 

.182 

.264 

.138 

.265 

.116 

.269 

20000 

.270 

.410 

.180 

.412 

.133 

.420 

50000 

.524 

.721 

.327 

.723 

.193 

.736 
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Table  4.8.  Computed  powers  when  the  RC  model  holds. 


Table  size 

3x3 

5x5 

10x10 

a N 

n2 

n4 

E2 

"4 

n2 

"4 

5000 

.010 

.011 

.010 

.013 

.010 

.011 

10000 

.011 

.012 

.011 

.017 

.010 

.012 

.01  20000 

.011 

.013 

.012 

.024 

.010 

.014 

50000 

.014 

.017 

.017 

.051 

.011 

.020 

100000 

.018 

.025 

.026 

.105 

.011 

.031 

5000 

.051 

.052 

.052 

.060 

.050 

.053 

10000 

.052 

.054 

.054 

.070 

.050 

.056 

.05  20000 

.054 

.058 

.058 

.091 

.051 

.061 

50000 

.062 

.071 

.073 

.154 

.052 

.079 

100000 

.074 

.093 

.099 

.262 

.055 

.108 

5000 

.102 

.103 

.103 

.115 

.100 

.104 

10000 

.103 

.106 

.107 

.129 

.101 

.108 

.10  20000 

.107 

.112 

.114 

.159 

.102 

.117 

50000 

.118 

.131 

.136 

.243 

.104 

.142 

100000 

.137 

.162 

.175 

.374 

.108 

.183 

CHAPTER  FIVE 
ISOTONIC  ESTIMATION 


Thus  far  we  have  presented  methods  for  the  analysis  of  categorical 
data  with  ordinal  variables.  This  includes  modeling,  interpretation, 
and  estimation  of  the  model  parameters.  If  there  is  a natural  apparent 
order  that  applies  to  the  levels  of  a categorical  variable,  instead  of 
assigning  arbitrary  monotone  scores  to  the  levels  of  that  variable, 
then  we  might  want  to  estimate  the  scores.  However,  limitations  of 
some  of  those  models  arise  when  the  scores  for  the  ordered  categorical 
variables  must  be  estimated.  For  instance,  the  estimated  scores  in  the 
RC  model  (2.5)  need  not  be  monotonic,  when  in  fact  the  property  of 
monotonicity  is  desirable  if  one  can  assume  the  true  scores  to  be 
monotonic. 

In  this  chapter  we  will  propose  procedures  for  two-way  cross- 
classifications that  produce  monotonicity  of  the  estimated  parameter 
scores.  A review  of  isotonic  regression,  the  primary  tool  used  to 
obtain  the  estimates,  will  be  outlined  in  section  1.  Sections  2 and  3 
deal  with  proposed  methods  for  obtaining  estimates  of  the  parameters 
for  the  R model  and  RC  model  respectively,  subject  to  the  constraint 
that  the  scores  be  monotonic.  An  investigation  on  the  performance  of  the 
proposed  procedure  of  section  2,  through  a simulation  study,  will  be 
discussed  in  section  4. 
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5.1 . Isotonic  Regression 

The  basic  theory  of  isotonic  regression  will  be  outlined  as  in 
Barlow  et  al.  (1972).  The  focus  will  be  primarily  on  definitions  and 
results  which  serve  our  purposes.  Let  E = {x] ,x2> . . . ,xk>  be  any  finite 
set  with  the  simple  order  x^<x2<. . .<xk- 

5.1.1.  Definitions 

A real  valued  function  on  E is  isotonic  (nondecreasing)  if  x,  yeE 

and  x<y  imply  f(x)<f(y).  Let  g be  a given  function  on  E and  w a given 

positive  function  on  E.  An  isotonic  function  g*  on  E is  called  an 

isotonic  regression  of  g with  weights  w if  within  the  class  of  isotonic 

functions  f on  E,  g*  minimizes  the  sum  S(g(xi )-f (xi ) )2w(x. ). 

In  what  follows  we  will  interpret  the  above  definitions  graphically. 

j 

Suppose  we  plot  the  cumulative  sums  G . = I g(x.)w(x.)  against  the 

j J i=l  1 1 

cumulative  sums  W.  = £ w(x.)  for  j = l,2,...,k.  That  is,  plot  the 
J i=l 

points  P.  = (Wj.Gj),  j = 0,1,2 k (define  PQ  = 0,0))  in  the  Cartesian 

plane.  These  points  connected  by  straight  lines  constitute  the 

cumulative  sum  diagram  (CSD)  of  the  given  function  with  weights  w. 

The  slope  at  the  point  P.  of  the  CSD,  approaching  the  point  from  below, 

is  g(xj ) = (Gj"Gj_i )/(Wj_Wj_i )•  The  9raph  of  the  supremum  of  all  convex 

functions  whose  graphs  lie  below  the  CSD,  is  referred  to  as  the  greatest 

convex  minorant  (GCM).  Graphically,  the  GCM  is  a path  along  which  a taut 

string  lies  if  it  joins  Pq  and  P^  and  is  contrained  to  lie  below  the 
★ 

CSD.  Let  P.  denote  the  point,  with  abscissa  W . , which  is  the  vertical 
J J 

projection  of  P^  on  the  path  constituting  the  GCM.  Then  the  slope  at 

* * 
the  point  P.  of  the  GCM,  approaching  the  point  from  below,  is  g (x.). 

J 
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5.1.2.  Preliminary  Results 

We  will  quote  two  results  from  Barlow  et  al . (1972,  pp.  12  and  64), 
which  will  be  used  subsequently. 

Theorem  5.1 

If  E is  simply  ordered,  the  slope  g*  of  the  CGM  furnishes  the 

isotonic  regression  of  g.  Indeed,  if  f is  isotonic  on  E then 

2 (g (x ) --f (x ) ) w(x)^z(g(x)-g  (x))2w(x)  + Z(g  (x ) -f (x) )2w(x).  The  isotonic 
xx  x 

regression  is  unique. 

Theorem  5.2 

Let  y be  an  unknown  function  on  a finite  set  E,  known  to  be 

isotonic  on  E.  Let  w(x),  xeE,  be  a set  of  positive  weights.  Let  g be 

* 

an  estimate  of  y.  Let  g be  the  isotonic  regression  of  g with  weight  w. 
Then 

T(y(x)-g*(x))2w(x)<z(y(x)-g(x))2w(x). 

X X 

The  latter  result  states  that  if  an  unknown  function  is  known  to 
be  isotonic,  then  the  isotonic  regression  of  an  estimate  of  that 
function  will  provide  a better  estimate  in  terms  of  mean  squared  error. 
The  former  result  will  be  used  to  show  the  uniqueness  of  the  isotonic 
regression  estimates. 

In  the  subsequent  sections,  we  use  the  definitions  and  results  above 
to  propose  methods  that  produce  isotonic  estimates  of  association  para- 
meters in  models  for  two-way  cross-classifications. 

5.2.  R Model 

Consider  a contingency  table  with  r rows  and  c columns,  with 
sample  size  N,  and  assume  one  of  the  three  sampling  schemes  (full 
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multinomial,  product  multinomial,  or  independent  Poisson).  The  R model 
is 


The  odds-ratio  for  the  four  cells  in  rows  i,  i+1  and  columns  k,£(k<£,)  is 


Although  the  R model  is  appropriate  for  a table  with  only  one  ordered 
categorical  variable,  we  still  can  use  the  model  for  a table  in  which  both 
variables  are  ordered.  In  other  words,  the  association  parameters 
{8.j}  can  be  treated  as  scores  rather  than  row-effects.  However,  the 
estimated  score  parameters  need  not  be  monotonic.  In  this  section, 
our  aim  is  to  estimate  those  parameters  subject  to  the  constraint  that 
the  parameter-scores  be  monotonic.  That  is,  we  want  to  fit  the  R model 
to  the  table  subject  to  the  following  constraint: 


= Npi j = exp  (y+A^+Xj+(vj-v)Bi). 


= exP(v£-vk)(Bi+1-3i)  l^isr-l , l<k<£<c.  (5.1) 


(5.2) 


Recall  that  to  find  the  ML  estimates  (m.  .}  of  the  expected  cell 
counts  under  the  R model,  one  must  solve  the  set  of  equations 
(2.10a-2.10c),  or  equivalently  the  following  set  of  equations: 


1 <i  <r 


lij^c 


1 <i <r  . 


(5.3) 
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We  refer  to  ft.  as  the  observed  or  estimated  row  mean  i under  the  R model 
for  i = 1,2,. ...r.  It  will  be  shown  later  that  the  estimated  row 
means  are  ordered,  if  and  only  if  the  estimated  B-parameters  are.  Hence 
we  converted  our  original  problem  to  a problem  of  ordered  row  means. 

That  is,  originally  we  wanted  to  fit  the  R model  subject  to  the  con- 
straint (5.2).  This  constraint  is  equivalent  to  the  constraint 

A A A 

R-|^R2<..  .<Rp.  Thus,  we  will  transform  the  table  in  a manner  which  will 

be  described  next,  such  that  the  observed  row  means  are  ordered,  and 

then  we  fit  the  R model  to  that  table.  Now,  we  are  ready  to  give  a 

rule  for  obtaining  estimated  expected  counts  {m*.}  under  the  R model 

■ *3 

subject  to  the  constraint  that  the  estimated  row-means  {R.}  are  ordered. 

Compute  the  observed  row  means  {R^ },  1 <i <r . Starting  at  k = 1, 
if  any  consecutive  pair  Rk,Rk+]  are  such  that  Rk>Rk+r  form  the  following 
weighted  average 


fi*  _ n*  _ nkA  +nk+l,A+l 

k - Rk+1  ' - V+nk+1>+  ' 


(5.4) 


It  will  be  shown  later  that  R^  and  R^  are  the  row  means  of  a trans- 
formed table  with  the  same  entries  as  the  original  table  but  rows  k and 
k+1  are  such  that 


n 


kj 


k+ 


V + Vi,f 


and 


n,  , . . = — 
k+1 , j n 


nk+l,+ 
k+  + nk+l ,+ 


(nko+nk+l,j>> 


(5.5) 


* * 

Note  that  n^+  = n.+  for  i = k,  k+1  and  n+j  = n+^,  l^j^c.  We  also  note 
that  rows  k and  k+1  have  identical  conditional  distributions  on  the 


104 


column  variable  in  the  transformed  table.  Now  we  have  r row-means  of 
which  r-2  are  unaffected,  with  row-means  R^+-|  being  equal.  Now 
start  from  row  k-1 , and  proceed  exactly  as  before,  treating  Rk=Rk+1  as 
a single  number  with  weight  (nk++nk+1  ),  and  treating  the  transformed 

table  as  the  original  one.  Continue  in  this  manner  until  the  observed 
row  means  are  ordered.  Finally,  we  obtain  a table,  each  row  of  which 
is  a linear  combination  of  the  rows  of  the  original  table,  that  has  all 
the  row  means  ordered.  The  fit  of  the  R model  to  the  final  transformed 
table  will  produce  ordered  estimates  of  the  6-parameters. 

The  proposed  solution  is  a special  case  of  isotonic  regression 
with  weights  equal  to  the  row  sample  sizes.  The  method  used  to  obtain 
the  solution  is  referred  to  as  the  "Pool  Adjacent  Violators"  algorithm 
in  Barlow  et  al.  For  the  general  notation  given  earlier,  the  finite 
set  E is  the  set  of  integers  {l,2,...,r},  the  function  g is  the  set  of 

A 

row-means  (g(i)  = , l^i<r),  and  the  weights  w are  the  row  sample 

sizes  (w(i)  = n..+,  l<i<r).  Then  the  isotonic  regression  g (i.e.: 

9 (i)  = R-j > l-i«r)  is  the  isotonic  function  that  minimizes  the  sum 

n A 12 

Z (R--R-)  n.  , over  all  isotonic  functions  { Rt } . 
i=l  1 

We  now  give  an  example  to  illustrate  the  method  graphically.  The 

data  in  Table  5.1  are  quoted  from  Plackett  (1972,  p.  24).  They  were 

obtained  in  the  course  of  a pilot  inquiry  into  the  conditions  in  which 

school  children  do  their  homework.  Children  are  classified  according 

to  the  conditions  under  which  homework  was  carried  out,  and  according 

to  the  teacher's  rating  of  the  quality  of  that  homework. 
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Table  5.1.  Classification  of  1019  children  according  to  conditions 
under  which  homework  was  carried  out,  and  the  teacher's 
rating  of  the  quality  of  that  homework.  (Each  scale  is 
graded.  A is  the  highest  rating.) 


Homework 

conditions 

A 

Teacher's 

B 

rating 

C 

Total 

Row-means 

<v 

A 

141 

131 

36 

308 

-.341 

B 

67 

66 

14 

147 

-.360 

C 

114 

143 

38 

295 

-.258 

D 

79 

72 

28 

179 

-.285 

E 

39 

35 

16 

90 

-.256 

Total 

440 

447 

132 

1019 

For  illustrative  purposes,  we  assume  that  the  scores  for  the 
categories  of  teacher's  rating  (grading)  are  chosen  to  be  integers:  -2, 
-1,0, 1,2.  We  also  assume  that  the  scores  for  the  categories  of 
homework  conditions  are  unknown,  but  known  to  be  ordered.  We  focus 
first  on  the  graphic  interpretation  of  the  isotonic  estimation  pro- 
cedure. 

We  see  from  Table  5.1  that  the  first  two  row-means  are  out  of  order 

as  are  the  next  two  row-means.  The  corresponding  transformed  row-means 

are  reported  in  Table  5.2.  In  this  table,  all  the  necessary  computations 

for  the  plotting  of  the  CSD  and  GCM  (Figure  5.3),  are  also  reported. 

★ 

To  compute  g by  the  algorithm,  note  that  g(l)>g(2),  so  we  form  the 
weighted  average  g*(l)  = g*(2)=^^^^°j+|~^60M147)  = .,347. 

Correspondingly,  the  first  two  rows  of  Table  5.1  are  transformed 
according  to  (5.5).  The  first  two  row-means  are  below  g(3),  so  next 
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we  should  modify  the  row-means  3 and  4 as  above.  That  is, 


9 (3)  = g (4)  = 


258  )(295 )+(-., 285 )( 179), 


(295  + 179) 


=-.268.  Therefore,  we  finally 


obtain  g (1)  = g (2)<g  (3)  = g (4)<g(5),  so  all  the  row-means  are  in 

order.  The  corresponding  transformed  table  has  its  row-means  ordered. 

Graphically,  since  g(l)>g(2),  the  part  of  the  GCM  between 
* * 

p0(=p0)  and  p2(=p2)  is  a straight  line  segment.  Thus,  the  CSD  could 
be  altered  by  connecting  pQ  and  p2  by  a straight  line  segment  without 
affecting  the  GCM.  As  seen  in  Figure  5.3,  similar  remarks  hold  for  the 
case  g (3)>g (4 ) . 

We  will  reconsider  the  example  later  and  will  give  the  transformed 
table  as  well  as  the  isotonic  estimates  {$*}  (p.  141). 


Table  5.2.  Example  of  CSD  and  GCM. 


j 

w(j) 

"j 

g(j) 

G. 

J 

g*(j) 

* 

gj 

1 

308 

308 

-.341 

-105 

-.347 

-106.876 

2 

147 

455 

-.360 

-158 

-.347 

-157.885 

3 

295 

750 

-.258 

-234 

-.268 

-236.945 

4 

179 

929 

-.285 

-285 

-.268 

-284.917 

5 

90 

1019 

-.256 

-308 

-.256 

-307.457 

* J * 

Note: 

g,-  = ^ g (j)w(j)  j 

1 ) • • • ) 5 

J i=l 
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Figure  5.3.  Graph  of  CSD  and  GCM. 

Before  deriving  properties  of  the  proposed  solution,  we  indicate 
a more  intuitive  new  approach  that  is  equivalent.  The  previous 
algorithm  suggests  that  we  first  compute  the  row-means  { R^. } to  determine 
their  ordering.  Since  the  { } are  in  the  same  order  as  the  estimated 
parameters  {f^-}  are,  it  is  then  equivalent  to  fit  the  R model  to  the 
rxc  table  to  obtain  the  ordering  of  the  (or  the  {ft.}).  Now  assume 
that  without  loss  of  generality  only  the  first  two  estimated  parameters 
are  such  that  B-| >&£  (°r  >^2 ^ * The  previous  method  suggests  that  we 

adjust  the  first  two  row-means  according  to  (5.4),  so  that  R*  = R*. 
Correspondingly  the  first  two  row-means  of  the  table  are  transformed 
according  to  (5.5).  The  fit  of  the  R model  to  the  transformed  table 
will  produce  B-|  = since  R-|  = R^.  Thus,  adjusting  the  first  two 
row-means  is  equivalent  to  fitting  the  independence  model  to  the  2xc 
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table  formed  from  rows  1 and  2,  since  el2  = 1,  l<j,j'<c  when  the  R 

J J 

model  is  applied  to  the  transformed  table.  If  the  adjustment  of  only 
the  first  two  row-means  yields  the  desired  ordering,  then  we  obtain 
(see  property  9)  that 

G2(IS0)  = G2(I)  + G2(Rco1),  (5.6) 

2 r c 

where  G (ISO)  = 2 E E n..  log(n../m..)  is  the  isotonic  goodness-of- 

i=l  j=l  1J  1J 

o 2 C j 

fit  statistic,  G (I)  = 2 E E n..  log(n../m!.)  is  the  statistic  for 

i=l  j=l  1J  1J  1J 


testing  independence  in  the  2xc  subtable  formed  from  rows  1 and  2,  and 
2 r-1  c 

G ^Rcol^  = ^ 2 2 nii  log(n. ./m. .)  is  the  goodness-of-fit  statistic 

i=l  j=l  J ,J  IJ 

for  fitting  the  R model  to  (r-l)xc  subtable  found  by  combining  rows  1 and  2. 

More  generally,  the  new  method  for  fitting  the  constrained  R model 
can  be  stated  as  follows: 

(i)  Compute  the  row-means  {R..}  and  adjust  them  according  to  (5.4). 
Suppose  we  obtain: 


R^R2^...<Rk 


/N* 

• • • = +■;  =•  • = 
K1  nl  k£ 


* = Rk  +i“’ 

£ £ 


.^R 


(ii)  Fit  the  independence  model  to  £ subtables,  each  formed  by  the 
sets  of  rows  for  lsu^£.  Let  G^(I  ) denote  the  goodness- 
of-fit  statistic,  for  testing  independence  in  subtable  u,  lSu^£. 


(iii)  Combine  in  a single  row  each  set  of  rows  for  which  the  in- 


dependence model  was  fitted.  The  resulting  collapsed  table  will  have 

ordered  row  means.  Let  G2(Rcq1)  denote  the  goodness-of-fit  statistic 

when  the  R model  is  applied  to  the  (r  - E i )xc  subtable.  By 

u=l  u 

performing  the  composite  fit  we  obtain  (property  9)  that 
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£ 

G2(IS0)  = G2(R  .)  + Z G2(I  ). 

C01  u=l  u 

Before  fitting  the  constrained  R model,  we  might  want  to  check  that 
lack  of  monotonicity  in  the  sample  is  not  that  severe.  For  instance,  if 
3-|>B2 (R-|>R2 ^ stron9ly  contradicts  Hq:  B-j  = b2,  then  we  will  interpret 

A 

the  { 3-j } as  row-effects  rather  than  scores.  Goodman  (1981)  discusses  how 
to  obtain  simple  tests  for  differences  between  association  parameters  in 


the  RC  model,  and  we  apply  the  method  to  the  R model.  Consider  an  rxc 
population  table  and  the  hypothesis  Hq:  61  = 82  in  the  R model.  Hence, 
the  R model  holds  and  Hq  is  true  if  and  only  if  both  of  the  following 
independent  conditions  hold 

(a)  There  is  statistical  independence  between  the  rows  and  columns 

in  the  2xc  subtable  formed  from  rows  1 and  2.  This  follows  from  the 
fact  that  e^2  = 1,  l<M<c. 

(b)  The  R model  fits  perfectly  the  (r-l)xc  subtable  found  by  collapsing 
the  first  two  rows  (see  property  1). 

2 

Let  G (I)  denote  the  chi-squared  statistic  for  testing  independence  in 


2 2 

the  2xc  table,  and  let  G (R)  and  G (R  .j)  denote  the  corresponding 
statistics  when  the  R model  is  applied  to  the  original  and  collapsed 
tables  respectively.  To  test  that  both  (a)  and  (b)  hold  simultaneously, 
we  can  use  G2(I)  + G2(Rcq1)  based  on  df  = (c-1)  + (r-2)(c-2). 

Under  the  assumption  that  the  R model  holds  in  the  rxc  table,  the 
statistic  G2(I)  + G2(Rcq1)  - G2(R)  based  on  df=[ (c-1 ) + (r-2) (c-2)] 

■ [ ( r-l ) (c-2 ) ] = 1 tests  Hq:  B-|  = B2-  The  method  presented  above 
can  be  extended  to  the  case  where  we  wish  to  test  the  hypothesis 


no 


that  a given  subset  of  the  parameters  3 ^ , l^iSr,  are  equal  to  each 
other. 

5.2.1.  Properties 

We  first  derive  two  properties  of  the  R model  and  by  symmetry  the 
C model  that  will  be  used  subsequently.  Then,  properties  of  the 
proposed  solution  will  be  derived. 


Property  1 

Suppose  we  want  to  fit  the  R model  to  an  rxc  cross-classification, 
and  suppose  that  the  first  two  rows  are  similar  in  the  sense  that 
one  is  a constant  multiplied  by  the  other  (i.e.,  n2j.  = an^,  l<j<c). 
Then,  the  fit  of  the  R model  to  the  rxc  table  is  equivalent  to  the 
fit  of  the  R model  to  the  (r-l)xc  table  found  by  collapsing  the 
first  two  rows,  in  the  sense  that 


(l+ajm^j  = (l+a_1)m2j. 


m 


i+l,j 


2£i2r-l 


where  the  {m9j}  denote  the  estimated  expected  counts  for  the 
collapsed  table. 

We  will  prove  the  claim  by  induction.  Consider  one  of  the  pro- 
cedures for  finding  the  ML  estimates,  say  the  iterative  scaling  method 

(Darroch  & Ratcliff,  1972).  Let  be  the  approximation  for  m. . 

th  ^ ^ J 

at  the  t stage.  Then  a single  cycle  has  3 steps 


m 


(t+2)  (t+1),  , (t+1 ) x 

mij  = V (niA-+  > 


m (t+3)  _ m (t+2) 
mij  'mij 


V** 

v . 
J 

zk(1'vk)nik 

„ * (t+2) 

^vkmik 

. . 

k * 

1-v. 


k = 0,1,2,. . . 


where  the  {vk>  are  a linear  rescaling  of  the  column  scores  {v^}  that 
. * c 

satisfy  0<vk<l . Let  {n^  l<i^r-l,  l^jSc}  denote  the  observed  cell 
counts  of  the  collapsed  table.  Note  that  n!j\  = (l+a^-jj  = (l+a_1)n2j, 
and  n^  = n^+-j  . for  2<i<r-l,  l^j<c.  Also,  the  collapsing  of  the  table 
does  not  affect  the  row  marginals  (i.e.,  n£.  = n .,  l<j<c). 

J 1 J 

To  fit  the  R model  to  the  collapsed  table,  we  start  with  m^.0^  = 1 

® J 

for  all  i,j.  Then  for  cycle  one  we  obtain 

il11  * "V(r-1) 


c (2 ) 

mij 


ni+n+j/N 


„c  . c 
c(3)  = Wj 

Vj 

jKKk 

mi  j N 

, 

N-'Ed-v^n^k 

1-v. 


(5.7) 


We  note  that  at  the  end  of  cycle  one  we  have 


c(3) 


- mi  ' 


i+1  »j 


2^i<r-l  and 


mf(3>  = 

(l+a)n1+n+. 

(HalEv*",. 

V . 
J 

(Ha)E(l-v*)nij 

lj 

N 

(Ha)Ev*n]+n+j 

• 

(l+a)EO-v*)n1+n+J 

= (l+a)m|^  = (l+a_1)m^ 


★ 

1-v. 

J 


(5.8) 


Assuming  that  (5.7)  holds  at  cycle  t,  we  prove  that  (5.7)  holds  at 
cycle  (t+1).  That  is,  at  cycle  t we  have 
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mljt}  " 0+a)mj^  = (l+a_1)m^  and  m^j^  = . for  2^i<r-l  and 


Then  at  cycle  (t+1)  we  obtain 


<St+,)  ■ . -&11  for  2£1£r-, , 


i+1 » j 


and 


'Ij 

we  have 


mljt+1)  = (^+a)m|j } (n+j/m|j  ^ ) = (l+a)m1(*+1)  = (1+a"1  )m^+1 } . Also, 


Ij 


2j 


mc!t+2> 

U 


™iJ(t+1)<">iit+1)> 


■ ■ *£!]  for 


2^i^r-l,  and 

»fJt+2)  - (Ho)"i1(‘+1)0+a)n..X(Ha)m.(‘+1))  = (Ha)»f?t2) 


m. 


] j V ' Ij  v U,"l+/u'  “;'"1+  ' 

= (1+a  ^ ' * Finally,  we  obtain 


m 


c(t+3)  _ c(t+2) 


ij 


= m.  : 
ij 


= m(t+2) 
mi+l , j 


* 

vj 


2(,-vk>4 


1-v. 


V . 
J 

2(1-vk)nik 

. , 

* 

’-vj 


= for  2<i^r-l,  and 
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m, 


c ( t+3 ) 

1j 


(l+a)m{j+2) 

0+a)Ev*nlk 

Vj 

(l+a)E(l-v*)nlk 

k 

<1+a>fkralk+2) 

h 4 

(l+a)E(l-v*)mjk+2^ 

k 

(l+a)m|j+3  ^ 

= (l+a_1)m^+3^  Q.E.D. 

Thus,  in  general,  when  we  fit  the  R model  to  an  rxc  cross-classification, 
we  can  collapse  over  the  rows  that  are  similar  in  the  sense  described 
above.  This  collapsing  will  not  affect  the  estimated  expected  counts, 
so  the  distances  between  the  estimated  B-parameters  are  preserved  in  the 
R model.  That  is,  if  we  let  the  {0^  and  {09}  denote  the  estimated 
association  parameters  respectively  for  the  original  and  collapsed  table, 
then  we  obtain 


flisi+1 

0k,£ 


A a 

mi  kmi+l  ,l  = Vl,kmU  _ 

/v  /S  /sQ 

mi£mi+l  ,k  mi-l,£mik 


= exp((v£-vk)(Bi+1-6i)) 


= exp((v£-vk)(§5-Bj_1)  2<i<r-l,  l<k<£<c; 


A A A A 

m,  ,m0n  am, . m. 


l o "h  « emu  . in,  . 

Ski  " ' 1 - exp[(B  -e  )(v  -v  )]  iskswe. 

mum2k  amumlk 


Hence  we  have:  §1  = §2  and  0.+1  - Bi  = 0*?  - 0^  1 2<i<r-l 


Property  2 

Suppose  we  want  to  fit  the  R model  to  an  rxc  table  for  which  the 
first  two  columns  are  similar  in  the  sense  that  one  is  a constant 
multiplied  by  the  other  (n^  = bn^,  l^i^r).  Then  assigning  the 
same  score  to  the  first  two  columns  and  fitting  the  R model  to  the 
rxc  table  is  equivalent  to  fitting  the  R model  to  the  rx(c-l) 


subtable  found  by  collapsing  the  first  two  columns,  in  the  sense 
that 
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ntjl  = (l+b)mil  = (l+b_1)mi2 


ai  j = ®u+i  2sJ«-i  • 

As  before,  the  claim  will  be  proven  by  induction  through  the 
iterative  scaling  method.  We  first  note  that  n^  = n^  + ni2  = (l+bjn^, 
and  ni j = ^ j+-j  2<j<c-l.  Also  the  row  marginals  are  unaffected  by 
collapsing  the  table  (i.e.,  n^+  = n^+,  l^i^r),  whereas  for  the  column 
marginals  we  have  n+]  = (l+b)_1n^,  and  n+  ^+1=  n^  for  2^j<c-l. 

At  cycle  one  of  the  iterative  process  we  obtain  (5.6)  and  at  the 
end  of  the  cycle  we  obtain 


-i?’  ■ 


■ * ^ 

Evknik 

vi 

EO-v*)n.k 

N',£vkni+"+k 

H-1zO-v*)n1+n+k 

1-V, 


Note 


that  Ivknn  = vi(nii+ni2> + k£2Vn  * Vii  + k22Vn  ' £Vn 


and  similarly  Z(l-vk)nlk  = E(l-vk)n^k,  and  thus  we  have 
k k 


* c 1 

Evk"1k 

V1 

£(1-vk)nik 

N'XniXk 

N-'jO-vk)n9+n^k 

”ii3)  ■ <1+b>'\% 


= (l+b)_1mfjy)  and  mj3^=  (1+b-1  )"1mj|3^ 
and  for  2Sj^c-l. 


1-v, 
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We  now  assume  the  claim  is  true  at  cycle  t,  and  then  prove  the  claim 
still  holds  at  cycle  t + 1.  That  is,  at  the  end  of  cycle  t we  have 

i-^  = (l+b)‘1m?jt\  and  m!^  = m^,  2SjSc-l - At  cycle  t + 1 we 
have  the  following: 


m 


m 


ilt+1)  = = (1+b)'1m-1(t)(n5+/mjjt))  = (l+b)"1n{{t+1). 


il  vi+/i+ 


11  v"i+/mi  + 


‘il 


and  m?2+1)  = (l+b-1  )'1m51(t+1 ) 


il 


} = milil(ni+/mH)  = mljt)<nl+/n^it>>  = raij(t)  ^c-l. 

This  follows  from  the  fact  that  = Zm^  = + m^  + Z m!^  = 

i+  1J  il  12  . >2  ij 

mil^^  + = mi+^‘  ^ seconb  sbeP  °f  the  cycle  we  obtain 


j>2 


">ilt+2)  = mi ] +1  ’ <n+l /mil+1  ’ )*  (Hb)‘1m^t+1)((Hb)-1n^/(l+b)-,m'{t+1))= 
and  m{*$  = 

(n+j/m+j^  ^ ) = Finally,  at  the  end  of  the  cycle  we  note  that 

£v{r>  s'r’) + * 

k>  i 

• and  simi,ariy  ^i-v><r2) = i'i'vk'mikt+2*' 


k>2vik  - 


Thus 


,(t+3)  _ 
"il 


(l+b)-1™?'4*2) 


(l+b)  and 


f * c 1 

zvk"ik 

V1 

fE(1-vk)nik  ] 

Zv*mc't+2> 
k lk 

1-v, 


m.(t+3) 

ray+i 


= mc<t+3> 


IJ 


=j=c-l 


Q.E.D. 
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More  generally,  suppose  we  have  an  rxc  cross-classification  which  con- 
tains sets  of  columns  which  are  similar.  The  property  states  that 
fitting  the  R model  to  that  table,  provided  the  same  score  is  assigned 
to  the  set  of  columns  that  are  similar,  is  equivalent  to  fitting  the 
R model  to  the  subtable  found  by  collapsing  the  columns  that  are  similar. 

Property  3 

If  an  rxc  cross-classification  has  only  two  distinct  columns,  then 
the  R model  fits  the  table  perfectly,  provided  that  only  two  dis- 
tinct scores  are  assigned  to  the  columns. 

In  this  situation  we  have  n^.  = a^n^ , l<j<s,  and  n^.  = b.^  s+1 


s+l<j<c. 

From  property  2 

above,  we  deduce  that 

( (a./  l a.  )m9 

ISj^s 

A 

m.  . 

_ ) J k<s  k 11 

i J 

l(Vk^sbk>*?2 

s<j<c . 

(5.9) 

Since  df  = (r-l)(c-2)  for  the  R model,  the  R model  is  saturated  when 
fitted  to  the  rx2  table  found  by  collapsing  the  first  s columns  and  the  last 
(c-s-1 ) columns,  i.e.,  m^.  = n|y  l^i^r,  j = 1,2.  From  equation  (5.9) 
we  conclude  that 


m.  . = 
ij 


((y,iakKi = ajnii = nij  isJss 


k<s 

vVJ/kKz = Vi.s+i = nij  s«£c- 


Hence  m^.  = n.^  for  all  i,j. 

Property  4 

The  final  set  of  transformed  row  means  { R_. } (consequently  the 
final  transformed  table)  through  equation  (5.4)  is  unique,  in  the 
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sense  that  we  obtain  the  same  result  in  the  algorithm  starting 
either  from  the  first  row-mean  or  the  last  one. 

This  property  follows  from  theorem  5.1,  since  the  {R*}  are  the 
isotonic  regression  of  the  { R^. }. 

Property  5 

As  mentioned  earlier,  transformation  (5.4)  is  equivalent  to  taking 

4>h  4-  U 

a linear  combination  of  the  k and  (k+1)  rows  and  then  computing 
the  row-means  of  either  row.  More  generally,  suppose  we  have  the 
following  situation 


A* 


R^R2^..  .5Rk 


* /\* 

Rk  +i  ^-^Rk 
K1 +11  k£ 


'"k  * 

• • • = Rk  +i  “* * *=Rr 

W r 


lSk1Sk,+1,S....SkJl  Sk^i.Sr 


Then  we  know  that 

k +i  k +i 

mm  mm 

Rk  = \ +i  = Z nx+Rx^  2 nx+’  1=m=5'» 

m x=km 


(5.10) 


that  is,  the  transformed  row-means  are  weighted  averages  of  the  original 
ones.  Now,  we  transform  the  corresponding  cells  of  the  original  table 
as  follows: 


k +i 
m m 


k +i 
m m 


/ I 
x=k 


m 


x+ 


"xj  for  km^”km+im*  1a"SK-  a"d  lsj«-  t5-11) 
m 


We  note  that  the  row  marginals  remain  invariant  under  (5.10),  since 


V ■ £ ■ V 


k +i 
m m 


j=l  x=k, 


x+ 


m 


k +i 
m m 

^k  V'V 

m 


fk  +i 
m m 

c 

Z In  . 
x=km  j=lXJ 


k +i 
m m 

'x*  "X+  = V. 

m 
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Now  the  row-means  of  the  rxc  table  of  which  some  cells  are  trans- 
formed according  to  (5.10)  will  have  its  row-means  exactly  as  the  ones 
in  (5.9).  Thus, 


Zv.n  ./n  , 

J yj  y+ 


£Vj(V 

?»Aj]/lBx+ 

X J 

fx+VZV  = Rk  = •••  = Rk+im  for  km"y^km+1m  and 
x m mm 


l=m=Jt-l . 

Hence,  in  general,  the  proposed  solution  is  equivalent  to  forming 
a new  table  whose  entries  are  linear  combinations  of  the  original  ones, 
and  whose  observed  row-means  are  ordered. 


Property  6 

The  column  marginal  totals  remain  invariant  under  the  proposed 
method. 

The  column  marginal  j of  the  final  transformed  table  is 

* 

n,.  = Z Zn  . + Z n . where  A is  the  set  of  integers  corresponding 
J m x XJ  xeA  XJ 

to  the  row  numbers  that  are  unaffected.  Now  Z Zn*. 


m x 


xj 


= 2 Z[n  (Zn  .)/Zn  ] = Z[(Zn  )(Zn xi)/(Zn.)]  = Z Zny..  Thus 

m x x XJ  x m x x XJ  x m x XJ 


Property  7 

When  the  R model  is  fitted,  the  estimated  (or  observed)  row-means 
of  an  rxc  table  have  the  same  order  as  the  estimated  B-parameters. 


Now  for  l^i^r-1,  we  have: 
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^ Vvk)mi  kVl  VV'V'Vl  ,*>° 


^/\-\)(aikai+urSiA+i,k)^/V\)Sik™ui,!1(§u1+1-1)>0 
«8kJ1+1>l.  lsk<«<c«-e1<Bi+r 


Property  8 

We  will  derive  the  asymptotic  distribution  of  6*  = (£*, . . . ,B*)T,  the 
estimate  of  produced  by  the  proposed  method. 

Let  _§  = (lp $2*  ’ * * *!r^  denote  the  ML  estimate  of  jB»  produced  by  the 

fit  of  the  R model  without  constraints.  We  refer  to  _§  and  §,  respectively, 
as  the  restricted  and  unrestricted  estimates  of  J3.  We  will  consider 
three  cases  separately. 

Case  1 

Let  us  assume  that  the  true-parameters  are  strictly  ordered  (i.e., 

. /S*  A 

^1 <^2< • * * <Br ) - Then  £ and  j3  will  have  the  same  asymptotic  distribution. 

Set  0.  = B,- , i -B.>0,  15i^r-l  and  choose  e such  that  0<e<6  = min  0.. 

1 =1=1 — 1 


Then  we  have 

P(®1+T®i>e»  1=i=r_1)  = P(Bi+1-Bi-6i>e-0i,  Ui^r-1) 

* P(|Bi+1-6i-ei |^0i-e,  l<i<r-l) 

= P(l  (3i+1-Bi+1  )-(Bi-6i ) l^e-Ei , l<i<r-l ). 
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Since  g is  a consistent  estimator  of  J5,  then 

p(|(Bi+1-8i+1)-(Bi-ei)|^ei-e,Ui^r-1)^P(|Bi-ei|^(ei-e)/2,1<i<r) * 1 

as  N 0°  . 

Thus  P(6i<62<- • *<6r)^P(B^+-| l£i£r-l)  ► 1 as  N + ».  That  is, 

VS>0  3NsjSuch  that  for  we  have  Pd^cB^. . . <3^ )>1  -6 . Using  property 

6,  we  conclude  that  for  N>N. 

0 

P(^1<P2<* • = P(8i<62<* •*<Bp)=l  - 6. 

Thus  for  a large  sample  of  size  NSNg,  the  row-means  will  be  ordered 
with  probability  at  least  1-6.  Hence,  as  N-*»  the  probability  converges 
to  1 that  the  original  table  does  not  need  to  be  transformed.  There- 
fore, both  the  restricted  and  unrestricted  estimates  are  identical. 

That  is,  V6>0,  3 N.  such  that  for  N>N„;  we  have 
o 6 

P(v^|l*-l|<e)>P(^B*=v#)=P(B*=B)=P(31<62  < — <3r)>l -6 
so  that  /N(|-l  ) — ^ 0. 

Since  v^CB-j?)  — N(0^,^]),  where  $ denotes  the  asymptotic  variance- 
covariance  matrix  of  the  ML  estimate  _B  (see  appendix  1),  and 
/N(j3-B  ) *■  0,  we  conclude  from  Slutzky's  theorem  that 

v¥(1*-B)  N(0,$). 

Thus  the  asymptotic  distributions  of  both  estimates  are  identical. 

Note  that  for  this  case  the  goodness-of-fit  statistic  for  the  con- 

strained  R model,  G (ISO),  will  have  the  same  asymptotic  distribution  as 

2 

the  goodness-of-fit  statistic  G (R)  for  the  unconstrained  R model. 

2 ? 

Both  statistics  G (ISO)  and  G (R)  will  have  an  asymptotic  chi-square 
distribution  based  on  df  = (r-l)(c-2). 


121 


Case  2 

For  this  case,  let  us  assume  that  p(2^p<r)  of  the  association 
parameters  are  out  of  order,  with  no  equality  between  them.  Note  for 
the  case  p = r,  we  will  have  only  one  distinct  row  in  the  transformed 
table,  and  hence  the  fit  of  the  R model  to  that  table  is  meaningless. 

For  simplicity,  consider  the  case  where  only  the  first  two  parameters 
are  out  of  order;  that  is  B-j >B2  but  B2<B3<. . .<B  . 

From  property  6,  we  know  that  only  the  first  two  population  row- 
means  would  be  out  of  order,  that  is  R,>R0  but  R0<R,<...<R  . If  we 
apply  the  proposed  algorithm  to  the  rxc  cross-classification  for  which 
the  first  two  row-means  are  out  of  order,  we  would  obtain  a population 
transformed  rxc  table  for  which  the  row-means  are  in  the  following 

k ic 

order:  R]  = Rg...  = RS<RS+1<- • -<Rr  for  l<s<r.  That  is,  the  final 

transformed  table  will  have  the  first  s rows  distributed  identically. 

If  the  R model  holds  for  the  transformed  table,  we  know  from  property  6 
* * * 

that  $i  = $2 • • • = ^s<es+i<* • *<3^  and  from  property  1 we  can  collapse 
the  first  s rows  without  affecting  the  expected  cell  counts.  Now,  we 
have  an  (r-s+l)xc  subtable  with  row-means  strictly  ordered  as  in  case  1. 
Let  J3  = (Bs»3$+i ,. . . ,$r)  denote  the  vector  parameter  associated  with 
the  collapsed  table,  and  let  BC  be  its  ML  estimate.  Then  from  Appendix  1 
we  conclude  that  /N(BC-BC)  — ^ N(0,Zc),  where  g^+1-^=6*  B*  for  s<i<r. 
Following  the  argument  of  case  1,  we  obtain  the  following: 

V 5>0,  3 N,.  such  that  for  N^Nr 

o 6 

P(^>R2  and  R2<R3<. . ,<R  )*l-6,  and 
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V*  /V* 


* 


P ( f^1  -^2_ * ' * RS<IW'--<V  = 


AW  A 


• < • • • <3^ )=1 _5 • 


If  we  collapse  the  first  s rows  of  the  transformed  table  with  sample 

size  N we  obtain  P(BSC<BS+-|<* • .<6rc)^l -6.  That  is,  for  N large  enough, 

with  high  probability,  the  subtable  of  size  N found  by  sampling  from  the 

original  table  and  then  collapsing  it  is  the  same  as  the  one  found  by 

sampling  from  the  collapsed  population  table  ( i . e . , P(B*c=j3c)>l -6  for 

N>N  ,). 

0 

Now  using  property  1 of  the  R model  we  obtain  that 
P(3i+1 -§*=£*+! -S*C,  s<i<r)>l-6  for  N>N6-  Hence 

f^[( (§i+1  “B.j ) - (Bi+1-Bi))  - - (b‘+1-b‘»]  0 for  s<i<r. 

Since  _B  is  asymptotically  normal,  then  from  Slutzky's  theorem  we  con- 
clude that  /tT[(6*+1-g*)  - (S*+1-S*)]  N(0,a*>i+1)  for  s<i<r,  where 

ai,i+l/N  is  the  asymPtotic  variance  of  B*+1  - g*.  We  will  show  that 

each  Bj  s^i<r,  is  asymptotically  normal.  To  do  so,  we  may  take 
r . 

A A "Jc 

.2  B-j  = 0 without  loss  of  generality,  so  that  vffir(Br-Br)  = 


JR 


/v*.  * * 

i=1 1 ^i+l"6i > - (ei+1_ej >> 


D . * * 2 

— * N(0,ar),  where  ar/r  N is  the 


asymptotic  variance  of  0 . Similarly  each  component  of  B will  have 


* 


an  asymptotic  normal  distribution.  Note  that  P(b1  =B2=- • -=BP H as  N-*». 


Hence  v¥(B*-B*)  N(0,t*), 


where  $ = 


* * 
Z11  Z12 

Z*T  Ec 
^12  L 


is  an  augmented  matrix  of  Z , the  asymptotic 
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covariance  matrix  of  _BC.  The  square  matrix  is  of  order  s-1, 

o^/N  is  the  asymptotic  variance  of  the  first  component  of  §c,  and  J 

is  an  (s-l)x(s-l)  matrix  of  ones.  The  matrix  Z1?  = 1 ,[ac  ,,...,aC  ], 

where  a^i+-|/N  is  the  asymptotic  covariance  of  6|+1  and  g?  for  s<i<r, 

and  ls_.j  is  an  (s-1)  vector  of  ones. 

Since  the  first  s rows  were  transformed,  the  isotonic  goodness-of- 

fit  statistic  G2(IS0)  = G2(Rcq1)  + G2(I),  where  G2 ( I ) is  the  usual  chi- 

square  statistic  for  testing  independence  in  the  subtable  formed  from  the 

first  s rows,  and  G (Rcq-j)  is  the  goodness-of-fit  statistic  for  fitting 

the  R model  to  the  subtable  found  by  collapsing  the  first  s rows, 
c c 2 

Since  3s<***<Br>  G (R^  ) will  have  an  asymptotic  chi-square  distribution 

o sc 

with  df  = (r-s+1 )(c-2).  The  statistic  G (I)  = 2 Z Z n..  log(n.  ./n. .n.  .) 

i=l  j=l  1J  1J  1+  +J 

= 2Nzf .j j 1 og (f i j/f i+f+j ) , where  f.^  = n^./N.  As  N+~,  f - £ 0 

where  f(£)  is  the  rc  vector  of  observed  (true)  cell  proportions.  Thus 
Zf -j j lo9(fij/fi+f+j)  D = Zp^-  log(pij./pi+p+j.),  as  N-*».  Since  the 

conditional  distribution  of  the  first  row  population  is  stochastically 
larger  than  the  conditional  distribution  of  the  second  row  population 
(3-j>62)’  we  know  that  D * 0.  Hence  G2(I)-*»  as  N+«>  with  probability  1, 
so  that  the  isotonic  goodness-of-fit  statistic  G2(IS0)-x»  with 
probability  1. 

More  generally,  suppose  that  p parameters  are  out  of  order.  If  we 
apply  the  proposed  algorithm  to  the  corresponding  rxc  cross-classification, 
we  would  get  a transformed  rxc  table  for  which  the  row  means  are  in  the 
following  order: 


R,  <R0<. . .<R.  =. 
1 2 k-| 


Kr 1 z r 


1,24 

As  for  the  case  of  only  two  parameters  being  out  of  order,  we 
would  collapse  the  rows  which  are  similar  to  obtain  a subtable  for 
which  the  population  row-means  are  strictly  ordered.  If  the  R model 
still  holds  for  that  subtable,  then  the  corresponding  vector  of 
association  parameters  will  have  an  asymptotic  normal  distribution; 
that  is,  /N(jsc-J3c)  — N(0,zc).  Hence  the  proposed  estimate  B*  will 
have  an  asymptotic  normal  distribution;  that  is,  /N(|*-B*)  — ^ N(jD,E*). 
As  for  the  case  p = 2,  the  asymptotic  covariance  E of  B is  an 
augmented  matrix  of  the  asymptotic  covariance  Zc  of  Bc. 

In  summary,  if  p parameters  are  out  of  order  for  an  rxc 
population  table,  we  apply  the  algorithm  to  the  table.  If  the  R model 
holds  for  the  transformed  population  table,  then  we  collapse  the  rows 
that  are  similar.  The  asymptotic  distribution  of  the  isotonic  estimate 
B.  is  obtained  through  the  asymptotic  distribution  of  Bc  for  the 
corresponding  collapsed  table. 

Case  3 

We  assume  here  that  p of  the  parameters  are  equal  (2^p<r). 

For  simplicity,  we  consider  the  case  where  only  the  first  two  parameters 
are  equal,  i.e.,  B-j  = ^2<^3<*  * *<&r*  For  this  case  1 = !*>  s0  that  the 
original  population  table  need  not  be  transformed.  In  what  follows  we 
will  present  a heuristic  proof  of  the  asymptotic  normality  of  the 
proposed  estimate  _B*. 

From  the  asymptotic  normality  of  B (see  Appendix  1)  we  obtain  that 
&l "-§2 ^ = (&] "£>1 )- ( )1  — N(0,a^),  where  o^/N  is  the 

asymptotic  variance  of  Thus  as  N -»■  °°  we  have  the  following: 
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P(8i<B3<--*<6r)  -*•  1 , i = 1,2  (same  argument  as  in  case  1),  and 
P ( B-j ^§2 ) = P(v/M(B-|-B2)  = 0)-^P( ) = i,  where  Z is  a normal  variate 
with  zero  mean.  Let  C]  = {61^62<63<* • -<Br}»  C 2 »{g  >g  ;B2<83<. . .4  }, 
and  let  C = be  the  complement  of  the  union  of  C-j  and  Let 

Al\|(l)  = {/N(6j-£)^z}  and  A^(zJ  = {/N(£-gJ<z}.  Then  we  have 

P{A„(i))  = P(C,  )P(A*(z)|C,  )+P(C2)P(A*(z.)|C2)+(1-P(C1  )-P(C2))P(A*(z.)|C). 


Let  and  4>^  (z^)  denote  the  asymptotic  normal  distribution  and 

density  functions  of  7n(£-&),  respectively.  Given  that  C-j  holds,  we 
know  that  | = £ so  that  P(AN(z)|C1)  = P(AN (^)  | C-, ) . For  a given  N 
large  enough,  we  know  that  P(C-j ) = { + o(l),  so  that 


“ $17(1)  + o(l),  where  the  derivative  of  is 
2c}) ^ (zj  over  the  region  C-j , 

'k0,  elsewhere. 


tIT(i)  -i 


Given  that  holds,  the  proposed  method  suggests  that  we  collapse 
the  first  2 rows  so  that  the  fit  of  the  R model  to  the  collapsed 
(r-l)xc  table  will  produce  estimates  such  that  •••<££  ]•  From 

case  1,  we  know  that  £c  = ($p§2> • • • »3^_]  )T  1S  asymptotically  normally 
distributed.  Let  and  4>2 ( z_)  denote  the  asymptotic  normal  distribu- 

tion and  density  functions  of  v¥(§a-^)  where  £a  = (3?,1C).  Hence,  given 

r * r zv 

that  C9  holds  and  X 6.  = E B?  = 0,  we  have  6 = (3a.  For  a given 

* i=l  1 i=l  1 “ “ 

N large  enough,  we  know  that  P(C2 ) = i + o(l),  so  that 
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p(AN(zJ|C2)  = + o ( 1 ) , where  the  derivative  of  $2T(z)  1S 

r 24>i  (z^)  over  the  region  C2 

<(>2t(I.)  ” | 

10,  elsewhere  . 

Therefore  it  follows  that,  as  N °°,  P(AN(z))  -*  i[$ITU)  + $2TU)]. 

Since  both  P(C-j)  and  P(C2)  tend  to  1 as  N goes  to  infinity, 

the  isotonic  goodness-of-fit  statistic  G2 ( ISO)  will  equal  the 
2 

statistic  G (R)  with  asymptotic  probability  | and  it  will  equal  the 
statistic  G2(I)  + G2(Rcol)  with  asymptotic  probability  £.  That  is, 
as  N oo 

P (G2 ( ISO ) = G2(R))  i and  P(G2(ISO)  = G2 ( I ) + G2(Rco]))  -4. 
Property  8 

If  the  population  row-means  are  ordered,  then  from  theorem 
5.2  we  know  that  the  restricted  estimate  row-means  provided  by  the 
proposed  algorithm  will  yield  estimates  with  smaller  mean  squared 
error  than  the  estimated  row-means  provided  by  the  standard 
method.  That  is,  l(Ri-Ri )2ni+<z(Ri-R*)2n.+  with  probability 
one. 


Property  9 

In  equation  (5.6)  it  was  stated  that  G2 ( ISO)  = G2(Rcq1 ) + G2 ( I ) , 
and  the  proof  proceeds  as  follows. 
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Assume  that  the  independence  model  was  fitted  to  the  subtable 
formed  by  the  first  s rows  ( i . e . , k-,  = 1,  i s-1  ).  In  this  situation 
the  transformed  observed  all  counts,  according  to  (5.5),  are  ex- 
pressed as  follows 

* 

n..  = n..,  i>s  and 

I J I J 

n*j  ‘ (ni+k^i nkj )/kE-1  nk+  * Vkj'  U1Ss- 

The  fit  of  the  independence  model  to  the  sxc  table,  produces 
estimated  counts  {m. .},  where 

• J 


-I 
m.  - 
U 


VVi 


s 

V 

^ nki 
k=l  KJ 

s 

Z n 
k=l 


■ "ij-  121ss- 


k+ 


The  fit  of  the  R model  to  the  data  set  {n..}  will  produce  estimated 

• J 

expected  counts  {m.^}.  Mow  if  we  collapse  the  first  s rows  of  the 
transformed  table,  we  know  from  property  1 of  the  R model  that  the 
estimated  expected  counts  {m'r.}  under  the  fitted  R model  to  the 

' J 

(r-s+l)xc  subtable  are 
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*ij=Vs-l,j  2^r's+1 


mC  - 


s 

1+  Z a. 
k=l  1 
k*i 


raij>  lsiss- 


Since  n.r.  = 

' J 


s 

1+  Z a, 
k-i 1 
k*i 


n..,  l<i<s,  then  we  have  n9 ./m9 . 
1 J i J 1 0 


n*j/iSij for 


l=i=s.  We  also  note  that  n..  = n . for  l>s  and  l<j<c.  Finally,  the 

' J I J 

three  statistics  in  equation  (5.6)  take  the  following  forms 
G2(IS0)  = 2Zn. . log(n../m*  ) = 2 z Z n..  log(n../m*  ) 

•j  'J  ij  j_i  j_i  lj  ij  ij 

c * 

+ 2 Z In..  log(n. ./m. .), 
i >s  j=l  1J  1J  1J 


s c 


GZ(I)  - 2 2 2 n,-j  log(n,./m!,)  -22  2 n.,  log(n. ./n. .)>  and 

1=1  j=l  J IJ  'J  1=1  j=l  'J  'J  ’J 


G <Rcol> 


2ln?j  ■ 2 2 

J 


r-s+1 


+ 2 2 nV  log(n=  /fig. ) 

.j_2  'J  iJ 


c 

= 2 Z 
j=l 


s 

Z n . 
i=l  1J 


* 


r-s+1 


"*og(nij/mij)  + 2 Z ns  . log(n..  ./m. . ) 


i=2 


ij 


ij  iJ 


s c 


* 


= 2,f1  ^,n1J  1og(nij/mij)  + 2 ,2,"ij  l09(nfi/ra(i), 


i>2 


n* 

i J*  ni j 


Thus  the  desired  result  is  obtained  by  writing 


G2(IS0)  - G2(Rco1)  - G2(I1)  = 2 


s c * s c 

Z Z n..  log(n../m z z n.. 
i=l  J=1  1J  1J  1J  i=l  j-1  1J 


j,"lj 


=0.  Q.E.D. 
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5.2.2.  Example 

For  the  data  in  Table  5.1,  we  fit  the  unconstrained  R model  and  the 

R model  subject  to  the  constraint  that  the  estimated  scores  be 

isotonic.  The  estimated  expected  counts,  the  transformed  cell  counts 

under  the  fitted  models,  the  goodness-of-fit  statistics  G2  and  the 

corresponding  estimated  association  parameters  are  displayed  in  Tables 

5. 3. a.  and  5. 3 . b . We  see  from  Table  5.3.b.  that  the  goodness-of-fit 
2 

statistics  G are  almost  identical  for  the  two  cases.  The  fit  of  the 
unconstrained  R model  suggests  that  the  odds  of  scoring  higher  in  the 
homework  decreases  in  the  order  of  the  homework  conditions  B,  A,  D,  C, 
and  E.  For  instance,  the  odds  of  getting  an  A rather  than  a B in  the 
homework  is  1.044  times  higher  for  students  who  carried  out  their 
homework  under  condition  B than  the  ones  under  condition  A.  If  one 
can  assume  that  the  row  categories  are  ordered  (with  A being  the  best 
condition,  E the  worst),  then  the  fit  of  the  R model  subject  to  the 
ordering  constraints  suggests  collapsing  the  first  two  row  categories 
and  the  next  two  row  categories.  That  is,  if  the  row  categories  A and 
B are  combined  as  well  as  the  row  categories  C and  D,  then  we  can 
interpret  the  sample  table  as  follows:  The  better  the  condition  under 

which  the  homework  was  carried  out,  the  higher  the  score  on  the  homework. 
For  instance,  the  odds  of  getting  an  A rather  than  a B in  the  homework, 
is  1.046  times  higher  for  students  in  categories  A or  B than  for  students 
in  category  C. 

To  illustrate  property  9,  consider  the  data  in  the  example.  Since 
the  first  two  row-means  are  transformed  as  well  as  the  third  and  fourth 
row-means,  we  fit  the  independence  model  to  the  subtable  formed  by  rows 
1,  2,  and  the  subtable  formed  from  rows  3 and  4.  We  obtain  the  goodness- 
of-fit  G ( I-j ) = .564  with  df  = 2 for  subtable  1 and  G2(I2)  = 3.117  with 
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df  = 2 for  subtable  2.  Collapsing  the  first  two  rows  as  well  as  rows 
3 and  4,  we  fit  the  R model  to  the  collapsed  3x3  table  and  we  obtain 
G2(Rcoi)  = 1-780  based  on  df  = 2.  Therefore  G2 ( ISO)  = 5.461  = G2 (Rcol ) 

+ G^(I^)  + G2 ( ) = 1.780  + 3.117  + .564.  When  the  R model  is  applied 
to  the  data  in  Table  5.3a,  the  likelihood  ratio  chi-square  G^(R)  = 

5.203  with  df  = 4.  To  test  the  hypothesis  H that  g.  = 60  and  g.  = gc, 
under  the  assumption  that  the  R model  holds  for  the  table,  we  can  use 
the  difference  between  the  two  statistics  G2 ( ISO)  - G2(R).  We  calculate 
the  difference  5.461  - 5.203  = .258  with  df  = 6 - 4 = 2.  Thus  we  do 
not  reject  the  hypothesis  H under  the  assumption  that  the  R model  holds 
for  Table  5. 3. a. 


5.3.  RC  Model 

In  this  section,  we  consider  the  case  where  both  ordered 
categories  in  the  rxc  cross-classification  have  unknown  scores.  The 
RC  model  is  naturally  designed  for  this  case,  since  it  has  two  sets 
of  score  parameters.  The  estimated  parameter  scores  produced  by  the 
fit  of  the  RC  model  need  not  be  monotonic.  However,  if  we  can  assume 
that  both  variables  have  unknown  yet  monotonic  scores,  then  we  can  fit 
the  RC  model  subject  to  the  constraint: 


A A 


A 


(5.12) 


If  we  let  S denote  the  subtable  which  contains  the  cells  in  rows 

{k, k+1 , . . . ,k+u}  and  columns  {£,£+1 , . . . ,£+v},  the  rows  and  columns  of 

S are  said  to  be  "quasi-independent"  if  m.  . = a.b.,  k^i<k+u, 

• J 1 J 
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£Sj££+v,  l<u^r-k,  l<v<c-£  (Bishop  et  al . , p.  178).  If  there  is 

statistical  independence  between  the  rows  and  columns  of  S, 

we  see  that  yk  = yk+1  = ...  = yk+u  and  VjL  = v^+1  = ...  = v£+v- 

To  test  quasi-independence,  denoted  by  Q,  in  subtable  S we  use  the 

chi-square  test  statistic  G2 (Q)  =2  In..  log(n . ./m9 . ) , and  the 

i,j  1J  1J 

summation  is  over  indices  of  rows  and  columns  in  S.  Expected 
frequency  estimates  {m9j } are 

™ij  = ni+(nkj+**'+nk+u,j)/(nk++--*+nk+u,+)  for  k^^k+u* 

™ij  = n+j(ni£+-“+ni,ii+v)/(n+£+'--+n+,ji+v)  for  £=j=*+v’  i<k  or  i>k+*> 


m. 


(n1+-Z 

X<£ 

X>£+V 


IX 


- I 

y<k 

y>k+u 


m . )/ 

yj" 


'k+u  £+v 

Z nxy 
x=k  y=£  y 


for  k<i<k+u  and 
£<j<£+v. 


Using  the  definition  and  the  remark  above,  we  propose  a method 
that  produces  isotonic  score  estimates.  For  ease  of  presentation,  we 
suppose  that  under  the  fitted  RC  model  the  estimated  row  scores 
Pk  anc*  ^+i  > as  well  as  the  estimated  column  scores  and  v^+^,are 

out  of  order.  Then  the  proposed  method  suggests  that 

(i)  We  fit  the  Q model  to  the  subtable  formed  from  the  union 
of  rows  k,k+l  and  columns  £,£+1,  and 
(ii)  We  combine  rows  k,k+l  and  columns  £ and  £+1,  to  obtain  a 
collapsed  table  for  which  the  RC  model  is  fitted. 

If  the  fit  of  the  RC  model  to  the  collapsed  table  produces  isotonic 
estimates,  the  proposed  estimates  are  such  that 
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y-,<y2^...<yk  = £k+1<. . .<£r>  and 


v1sv2<...<vl  = Vj+1<...<vc. 


Those  estimates  are  found  through  the  estimates  (pi),  (GC;  produced  by 

J 

the  fit  of  the  RC  model  to  the  collapsed  table,  using  the  relation 


A A A A A ^ A Q 

h+th  = ^i+i^i 


vj+rvj  = vj+rvj’ 


l£i<r-l , 


l^j^c-1  and  treating 


vWw,> 


as  a single  number. 

If  the  fit  of  the  RC  model  to  the  collapsed  table  does  not  produce 
isotonic  estimates,  then  we  combine  rows  and  columns  for  which  the 
corresponding  estimates  are  out  of  order,  and  repeat  the  procedure  described 
above.  Let  GRC( ISO)  = 2Zn^  logfn^/m^)  denote  the  goodness  of  fit  statis- 
tic. Then  it  can  be  shown  with  the  same  argument  used  to  prove  (5.6),  that 


grc(iso)  = g2(q)  + g2(rcco]), 

2 

where  G (Q)  is  the  statistic  for  testing  the  Q model  in  the  corresponding 
table,  and  G (RCcq1)  is  the  statistic  for  fitting  the  RC  model  to  the 
collapsed  table. 

As  for  the  R model,  we  can  check  whether  lack  of  monotonicity  in  the 
sample  is  not  that  severe.  For  instance,  if  yk>yk+1  and  v£>v£+1  strongly 
contradicts  Hq:  (yk  = yk+-|  and  v£  = v£+-|),  then  we  need  not  fit  the 
constrained  RC  model.  That  is,  it  makes  more  sense  to  interpret  the  scores 
as  effects  in  that  table. 

Goodman  (1981)  discusses  how  to  test  H . The  RC  model  perfectly  fits 
the  rxc  population  table  and  Hq  is  true  if  and  only  if  both  of  the 
following  conditions  hold: 
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(a)  There  is  statistical  independence  in  the  table  formed  from  the 
union  of  rows  k,k+l  and  columns  £,£+1. 

(b)  The  RC  model  holds  for  the  (r-l)(c-l)  collapsed  table  in  which 
rows  k and  k+1  as  well  as  columns  £ and  £+1  are  combined. 

Conditions  (a)  and  (b)  can  be  tested  by  G2(Q)  and  G2(RCcol ) respectively. 

To  test  both  conditions  (a)  and  (b)  simultaneously,  one  can  use  the  sum 
G2(Q)  + G2(RCcq1).  Hence  the  statistic  G2(Q)  + G2 (RCcq1 ) - G2(RC)  based 
on  df  = 2 tests  Hq.  The  method  presented  here  can  be  extended  to  the  case 
where  we  wish  to  test  HQ : (yk=Pk+1=. • .=Pk+u  and  v£=v£+1  = . . . =v£+y ) . 

5.3.1.  Properties 

We  first  derive  a property  of  the  RC  model  which  we  will  use  sub- 
sequently to  derive  the  following  properties. 

Property  1 

Consider  an  rxc  cross-classification  for  which  a set  of  rows  and/or 
a set  of  columns  are  similar.  Then  the  fit  of  the  RC  model  to  that 
table  is  equivalent  to  the  fit  of  the  RC  model  to  the  table  found  by 
combining  similar  rows  and  similar  columns. 

The  RC  model  can  be  fitted  to  the  table  by  fitting  alternatively 
the  R model  and  C model.  The  claim  stated  above  follows  from  properties 
1 and  2 of  the  R model  and  by  symmetry  those  of  the  C model . 

Property  2 

If  the  final  collapsed  table  for  which  the  RC  model  is  fitted  contains 
only  two  rows  (columns),  then  the  isotonic  goodness-of-fit  statistic 
equals  the  goodness-of-fit  statistic  for  testing  "quasi-independence"  or 
independence  between  the  rows  and  columns  of  the  corresponding  table. 
This  follows  from  the  fact  the  residual  degrees  of  freedom  for  the  RC 
model  df  = (r-2)(c-2)  = 0,  so  that  G2(RCco1)  = 0. 


134 


Property  3 

The  proposed  method  is  equivalent  to  finding  a transformed  table  such 

that  the  RC  model  produces  isotonic  estimates  when  applied  to  that  table. 
This  follows  from  property  1 and  the  fact  that  the  combined  rows  and 
columns  in  the  process  are  similar. 

Property  4 

^ ^ 

We  will  derive  the  asymptotic  distribution  of  _§  = (y-j  .....y^v-] , 

/v*  T T 

• • • >v  ) , the  estimate  of  the  vector  parameter  0 = (y,  ,...,y  ,v-,,...,v  ) , 

produced  by  the  algorithm  under  the  assumption  the  true  population 

scores  {y..}  and  {v j } are  strictly  ordered.  If  we  let  ^ = (y-j 

yr»v-| ....  ,vc)T  denote  the  ML  estimate  of  0 produced  by  the  fit  of  the 

unconstrained  RC  model,  then  ^ and  j0  will  have  the  same  asymptotic 

distribution. 

To  obtain  the  asymptotic  distribution  of  §*,  we  follow  the 
argument  of  case  1 in  section  2,  from  which  we  can  deduce  that 
P (y-j <y2< • • • <Pr ;v-| <v2< • • • <vc H-l  as  N-*».  That  is,  for  a large  sample 
of  size  N>Ng,  both  sets  of  estimated  scores  {y.}  and  {v.}  will  be 
ordered  with  probability  at  least  1-6.  Then  both  estimates  _§*  and 
_0  will  coincide  with  high  probability,  so  use  of  the  algorithm  is 
unnecessary.  Now,  as  N-*»  we  have  P(_§=0*)  = P(v^§  =/N0*)-*l,  so  that 
/N((j^  -_0)  - (0-j)))  — 0.  Since  § has  an  asymptotic  normal  distribution 
(see  Appendix  1),  from  Slutzky's  theorem  we  conclude  that 

/N((j3*-  0)  N(0,Z),  where  £ is  the  asymptotic  covariance  matrix  of  9. 

In  this  case  where  the  true  parameter  scores  are  strictly  ordered,  the 
isotonic  goodness-of-fit  statistic  will  have  an  asymptotic  chi-square 
distribution  with  df  = (r-2)(c-2)  as  the  regular  goodness-of-fit  statistic 
does  when  the  RC  model  is  applied  for  the  rxc  cross-classification. 
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5.3.2.  Example 

For  the  data  in  Table  5.4  on  dumping  severity  for  ulcer  operations, 
we  fit  both  the  RC  Model  and  the  RC  Model  subject  to  constraints  (5.12). 
The  parenthesized  values  are  the  estimated  expected  counts  under  the 
RC  Model  and  the  constrained  RC  Model,  respectively. 

Maximum  likelihood  estimation  of  the  RC  Model  gives  G^c  = 2.86  based 

on  df  = 2.  There  is  little  improvement  in  fit  over  the  uniform 

2 

association  model,  for  which  Gy  = 4.59  based  on  df  = 5.  The 
latter  model  gives  a good  fit,  and  suggests  that  dumping  is  more  severe 
when  more  of  the  stomach  is  removed.  However,  for  the  RC  Model,  we 
obtain  y-j  = .364,  y2  = -.613,  y^  = .404,  and  y^  = .573  for  estimated 
scores  on  operation  and  v1  = -.767,  v2  = .552,  and  v3  = .245  for  the 
estimated  scores  on  dumping  severity.  Since  the  estimated  scores  are 
nonmonotonic,  there  is  indication  of  nonmonotonic  associations,  in  the 
sense  that  local  associations  are  positive  in  some  locations  and 
negative  in  others.  For  instance,  the  odds  of  not  having  dumping 
rather  than  slight  dumping  severity  is  .72  times  as  high  as  for  patients 
who  had  operation  A than  for  patients  who  had  operation  B.  However, 
the  odds  of  having  slight  dumping  severity  instead  of  moderate  dumping 
severity  is  1.08  times  as  high  for  patients  who  had  operation  A than 
for  those  who  had  operation  B. 

To  use  the  ordinal  nature  of  the  categories  by  assuming  the  true 
row  and  column  scores  are  monotonic,  we  fit  the  constrained  RC  Model 
to  Table  5.4.  Since  the  estimated  row  scores  y^  and  y2  as  well  as 
the  estimated  column  scores  are  out  of  order,  we  fit  the  quasi- 
independence model  to  the  table  formed  from  rows  1,  2 and  columns 
2 and  3.  Then  we  fit  the  RC  Model  to  the  3x2  sample  table  by  combining 
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the  first  two  rows  and  the  last  two  columns.  We  obtain  that  the 
statistic  for  testing  quasi-independence  is  G (Q)=3. 03  with  df  = 4, 
and  G2(RCcq1)=  0.0  (property  1)  so  that  G^c( ISO)  = G2(Q)  + G2(RCcol)  = 

3.03  with  df  = 4.  We  also  obtain  ^ = P2  = --493,  y3  = .379,  and 

/v*  /v* 

P4  = .607  for  the  row  scores,  v-j  = -.816  and  v2  = v3  = -408  f0*''  the 
column  scores.  Since  the  estimated  scores  are  monotonic  nondecreasing, 
we  conclude  that  there  is  a tendency  for  dumping  to  be  more  severe 
when  more  of  the  stomach  is  removed.  However,  we  notice  that  the 
conditional  distributions  of  the  first  two  row  variables  are  identical 
as  well  as  the  conditional  distributions  of  the  last  two  column 
variables.  To  test  the  hypothesis  that  y1  = y2  and  v2  = v3  under  the 
assumption  that  the  RC  model  holds,  we  use  the  difference  between  the 
statistics  GRC( ISO)  and  GRC-  We  calculate  the  difference  3.03  - 2.86  = .17 
with  df  = 4 - 2 = 2.  Thus  we  cannot  reject  the  hypothesis  that  = y2 
and  v2  = v3  under  the  assumption  that  the  RC  model  holds.  Both  the 
constrained  and  unconstrained  RC  model  fit  the  data  quite  well. 

5.4.  Simulation  Study 

As  noted  earlier  (section  2,  property  7),  if  the  population  row- 
means  { R^ } are  truly  ordered,  then  the  proposed  method  for  the  R model 
would  produce  row-mean  estimates  {R^ } that  are  better  than  the  estimated 
row-means  (R^ } found  by  just  fitting  the  R model,  in  a mean  squared 
error  sense.  However,  we  could  not  derive  a similar  property  for  the 
estimated  scores  produced  by  the  two  different  procedures.  We  instead 
conducted  a simulation  study  to  investigate  the  performance  of  the 
estimated  parameters  { 3^ } , comparing  them  to  the  regular  estimates 
(Bi  }• 
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5.4.1.  Description 

In  this  simulation  study  we  consider  only  3x3  cross-classifications. 
For  all  the  tables,  each  row  is  an  independent  normal  variate  with  a 
certain  mean  and  with  unit  standard  deviation.  We  discretized  the 
normal  variate  X at  the  points  -.5  and  .5  as  follows:  the  cells  of 

a row,  say  i,  are  p^  = P(X<-.5),  pi2  = P(-.5<X<.5),  and  pi3  = P(X>.5). 

We  considered  two  different  shifts  A (.1  and  .5)  for  differences 
between  means  in  adjacent  rows.  The  first  shift  would  result  in  a 
small  spacing  of  the  3-parameters,  whereas  for  the  other  shift  the 
spacings  are  somewhat  larger.  For  each  shift  we  consider  two  cases. 

In  the  first  case  all  rows  have  different  means,  whereas  in  the  second 
case  the  last  two  rows  are  identically  distributed.  The  latter  case 
results  in  an  equality  between  the  last  two  3-parameters  in  the  R model. 
For  each  of  the  four  cross-classifications,  and  for  two  sample  sizes, 
we  fit  the  R model.  The  model  gives  a very  good  fit,  though  not  a 
perfect  fit.  Finally,  to  have  a population  table  for  which  the  R Model 
fits  perfectly,  we  instead  consider  just  the  tables  with  cell  entries 
being  estimated  expected  counts  under  the  fitted  R model.  That  is, 
for  all  the  3x3  cross-classifications  under  consideration,  the  R model 
fits  perfectly.  Furthermore,  all  association  parameters  are  known  to 
be  monotone  nondecreasing. 

The  purpose  of  this  simulation  is  to  estimate  n-  = E(B- -3- )2  and 
D-j  = E(3.j-3.j)  for  l<i<3.  To  do  so,  for  a given  sample  of  size  N (30 
and  150),  we  randomly  generated  K = 1000  3x3  contingency  tables,  where 
each  row  in  each  table  is  an  independent  multinomial  sample  of  size  N 
(10  and  50).  For  each  of  those  K tables  we  fit  the  R model,  then  the 
R model  subject  to  the  constraint  that  the  estimated  association 
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parameters  {6^ } be  monotonic.  In  simulation  £,  let  and  B*^ 

l^i^3,  denote  the  estimated  B-parameters  produced  by  their  respective 

★ 

method.  Then  the  estimates  of  and  n..  l<i^3,  would  be 


A 


l (B1Jt- Bt) 
£=1  u 1 


2 


,-l 


and  r\.  = K Z (8.-.-BJ  . 
£=1  * 


To  measure  the  performance  of  the  proposed  estimates  in  terms  of 
mean  squared  error,  we  look  at  the  ratio  of  the  mean  squared  error  of 
the  two  different  estimates  n*/n . for  1 <i<3  and  the  ratio 


3*3 
£ W £ r).. 
1=1  1 i=l  1 


Other  useful  statistics  have  been  reported  and  are  given 


next. 

(i)  The  estimated  expectations  {E(Bn- and  {E(B*)>, which  are  denoted 

★ K K 

by  {B.j}  and  fiL}  respectively,  where  B,  = K_1  z B-  and  'b*  = K_1  z 

£=1  U 1 £=1  U 

for  l^i<3. 

(ii)  The  approximation  of  the  } by  their  asymptotic  theoretical  values, 
using  the  asymptotic  normality  of  the  ML  estimate  JL  Since 

*^(_t-j5)  N(04),  one  could  approximate  n..  by  a2,l<i^3,  where  a?  is 

taken  from  the  diagonal  of  $,  i.e.,  n-  = E(B1--6i  )2=a2  Ui<3. 

(iii)  The  asymptotic  standard  deviations  of  the  estimates  (fj . } : 

Since  | is  asymptotically  normal,  then  E(Bi-3i )4 = 3a|  where  a?  is 
defined  as  above,  and  var[ (§i-Bi )2]  = 3a4  - a4  = 2a4,  l<i^3.  Hence 

K K 

we  have  <£  = var(n.)  = var[K_1  Z (B.  -B,)2]  = K"2  Z var(B, 0-B- )2=2a^/K, 

'i  £=1  1 £=l  1£  1 1 

and  a-  = y^/Ka2  1 <i <3 . 
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(iv)  The  sample  standard  deviations  of  both  estimates  {f^.}  and  {rj*} 

ys  * 

which  are  denoted  by  {o~  } and  {o^. } respectively,  where 


2/K  with  ni£  = (0i£-Bi)2  and 


with  ni£  = for  1 <i <3 . 


5.4.2.  Results 

We  will  present  the  results  of  this  simulation  study  in  two  separate 
cases.  We  consider  first  the  case  where  the  rows  have  different  dis- 
tributions, then  the  case  where  the  last  two  rows  are  identically 
distributed. 


For  the  two  shifts  A (.5  and  .1),  the  underlying  row  variables  for 
the  3x3  tables  are  independent  normally  distributed  with  unit  variance 
and  respective  means  -A,  0,  and  A.  For  this  case,  the  association 
parameters  {f^}  are  evenly  spaced.  The  spacing  increases  as  the  shift 
increases.  For  A = .5  we  obtain  g-|  = -g^  = -.576  and  g ^ - 0,  whereas 
for  A = .1  we  have  g^  = -g^  = -.104  and  = 0*  The  population  tables 
considered  in  this  study  are  reported  in  Table  5.5. 

From  Tables  5. 6. a. -5.6. b. , we  see  that  both  estimates  1$  and  g* 
provide  biased  estimates,  for  both  shifts,  of  the  first  and  third  com- 
ponent of  _g.  The  bias  decreases  as  the  sample  size  increases.  We  know 

that  the  {g£}  and  (g£)  for  1^£<K,  will  tend  to  {g}  in  probability  as 
★ 

N-*»,  so  that  _g  and  _g  will  tend  in  probability  to  g.  To  measure  the 


Case  1 
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accuracy  of  _g  and  jj  we  could  use  their  respective  sample  variance  ti. 
and  r\..  For  instance,  for  the  case  of  a small  shift  and  small  sample  size, 
we  obtain  = -.107,  g2  = *003,  and  g^  = *104  with  their  sample  standard  devia- 

tion {/n^lof  order  .380,  and  g^  = .208,  g^  = .218  with  sample  standard 

— * 

deviation  of  order  .300  and  g2  = -.010  with  sample  standard  deviation 
y/f)2  = .177.  In  this  simulation  study,  the  proposed  method  (to  a lesser 
degree  the  competitor)  underestimates  the  first  parameter  g-j , and  over- 
estimates the  third  parameter  g3*  In  this  case,  the  population  tables  have 
a special  form  in  the  sense  that  by  reversing  both  the  rows  and  columns 
they  remain  invariant.  Let  be  a given  sample  table,  and  let  be 
the  table  found  by  reversing  both  the  rows  and  columns  of  T(1  \ Assume  that 
T^1  produces  estimates  such  that  gj^>g^  and  g^<g^.  If  we  apply  the 
algorithm  to  7^1  we  obtain  estimates  such  that  gj^*=  g^*<g^*.  Now 
the  R model  applied  to  will  produce  estimates  such  that  g{2^<g^  and 
^22)>^32)’  where  g^=-g^\  and  g^^-gj1)'  The  algorithm 

applied  to  7^  will  produce  estimates  such  that  gj2^*<g22^*=g32)*,  where 
gj  =-63  ' and  g£  ‘ =-gj  ' . Hence  if  a set  of  tables  {T|  0 produce 
estimates  g2_.  and  g^,  then  the  corresponding  set  of  tables  {Tl  '}  will 
produce  estimates  -g^  and  -62-j  * Since  each  table  tP  ^ occurs  with  the 
same  probability  as  the  corresponding  table  t|2^,  the  distribution  of  fl, 
and  g2  are  symmetric  around  0.  Thus  both  methods  provide  an  unbiased 
estimate  of  the  middle  category  parameter  g^. 

In  terms  of  mean  squared  error  the  proposed  method  performs  better. 

The  performance  is  relatively  better  for  small  shifts  and  small  sample 
sizes.  It  is  already  known  from  property  7 of  the  R model  that  the  two 
methods  should  perform  equally  as  the  sample  size  gets  very  large.  Also 
in  this  case,  the  two  methods  should  perform  as  well  for  larger  shifts. 


A.  A* 

To  measure  the  preciseness  of  the  MSE  estimates  {q..}and  {r^}  we  use 

their  asymptotic  standard  deviations  {a~  }.  For  instance,  for  the  case 

nl 

of  a large  shift  we  obtain  rL  = .162,  ri0  = .103  with  o * = .005  for  a 

A /S* 

sample  size  N = 30,  and  we  have  n0  = .025,  n0  = -024  with  o * = .001 

C £.  T)^ 

for  the  larger  sample  size  N = 150.  For  small  N,  there  is  clear  evidence 
that  the  constrained  R model  produces  better  estimates  than  the  uncon- 
strained R model  in  mean  square  error  sense.  The  proposed  method  performs 
relatively  better  for  the  smaller  shift.  Finally,  we  remark  that  the 
sample  standard  deviations  {a^}  and  {a~*}  approximate  poorly  the 

asymptotic  standard  deviations  {cr'  },  for  small  sample  size.  The 

ni 

approximation  improves  as  N increases.  For  further  specific  comparisons 
of  the  two  methods  see  Tables  5. 6. a. -5. 6. b. 

Case  2 

For  the  two  different  shifts,  the  last  two  rows  have  the  same  dis- 
tribution. That  is,  the  first  row  is  distributed  as  in  case  1,  and 
the  last  two  rows  are  i.i.d  standard  normal.  For  this  case,  we  know 
from  property  6 of  the  R model  that  there  are  only  two  distinct 
B-parameters.  When  the  shift  A = .5,  we  obtain  B1  =-.384,  b2  = $3  = *192 
and  when  A = .1,  we  have  B-j  = -.068,  B2  = 83  = .034.  As  in  the  previous 
case,  the  estimates  produced  by  the  constrained  R model  yield  smaller 
mean  square  errors  than  the  estimates  produced  by  the  unconstrained  R 
model.  Actually,  we  see  from  Tables  5. 7. a. -5. 7. b.  that  the  performance 
of  the  proposed  estimates,  in  terms  of  mean  squared  error,  is  better 
for  this  case  than  the  previous  one,  particularly  for  the  middle  cate- 


gory parameter  B2* 
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In  summary,  we  would  expect  E(3*)<E(£-|)  and  E(£*)>E(B3) . In 
other  words,  in  simulation  £,,  suppose  we  fit  the  R model  to  the  sample 
3x3  table,  and  the  estimated  parameters  {§.,,}  are  out  of  order.  The 
proposed  solution  tends  to  lower  the  estimated  value  6^  whenever  the 
first  two  estimated  parameters  are  out  of  order.  It  also  tends  to 
increase  the  value  of  B3^  whenever  the  last  two  estimated  parameters 
are  out  of  order.  However,  for  the  middle  category  parameters  62>  the 
method  tends  to  increase  (decrease)  the  values  of  the  estimates  } 
whenever  the  first  (last)  two  estimated  parameters  are  out  of  order. 

We  also  would  expect  the  constrained  R model  to  produce  better 
estimates  than  the  unconstrained  R model  for  small  sample  sizes  and 
small  shifts.  To  examine  the  behavior  of  the  proposed  method  as  a 
function  of  sample  size,  shift,  and  nature  of  the  distributions  of  the 
row  variables,  we  computed  the  relative  sizes  of  MSE  (i.e., 

Zn^/Zn^ ) which  are  displayed  in  Table  5.8.  We  see  from  the  table  that 
the  relative  size  of  MSE  increases  as  the  size  of  the  sample  and/or  the 
shift  increase.  For  fixed  sample  and  fixed  shift,  the  performance  of 
the  proposed  estimates  is  better  in  the  second  case  than  in  the  first 


case. 
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Table  5. 3. a.  Cross-classification  of  1019  children  according  to 
conditions  under  which  homework  was  carried  out. 


Homework 

conditions  ABC 


A 

141 

(140.61) 

131 

(131.77) 

36 

(35.612) 

140.8 

(141.912) 

133.354 

(131.13) 

33.846 

(34.958) 

B 

67 

(69.028) 

66 

(61.943) 

14 

(16.029) 

67.2 

(67.730) 

63.646 

(62.585) 

16.154 

(16.685 

C 

114 

(118.948) 

143 

(133.103) 

38 

(42.949) 

120.116 

(120.851) 

133.808 

(132.338) 

41.076 

(41.811) 

D 

79 

(75.238) 

72.192 

(79.524) 

28 

(24.238) 

72.884 

(73.330) 

81.192 

(80.300) 

24.924 

(25.370) 

E 

39 

(36.174) 

35 

(40.652) 

16 

(13.174) 

39 

(36.176) 

35 

(40.647) 

16 

(13.174) 

Note:  Each 

cell  has  two  entries, 

the  observed  counts  anc 

i transformed 

counts.  The 

parenthesized  values 

are  the  estimated  expected  counts  for 

the  row-effects  model . 

Table  5 . 3 . b . 

Goodness- 

of-fit  test  statistics 

i and  estimated  association 

parameters  for  the 

R model . 

Model 

G2 

1 

2 3 

4 

5 

R model 

5.203 

A 

6 

-.087 

131  .090 

.033 

.095 

Constrained 

R model 

5.461 

B 

-.107 

107  .063 

.063 

.088 
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Table  5.4.  Cross-classification  of  duodenal  ulcer  patients  according 
to  operation  and  dumping  severity  (parenthesized  values  are 
estimated  expected  frequencies  corresponding  to  the  RC 
model  and  the  constrained  RC  model). 


Dumping 

severity 

Operation 

None 

Slight 

Moderate 

Total 

A 

61 

(60.28,  61.42) 

28 

(25.58,  24.84) 

7 

(10.12,  9.42) 

96 

B 

68 

(68.62,  67.08) 

23 

(25.10,  26.97) 

13 

(10.28,  10.01) 

104 

C 

58 

(57.60,  50.00) 

40 

(38.63,  37.90) 

12 

(13.77,  14.10) 

110 

D 

53 

(53.50,  53.00) 

38 

(39.68,  39.36) 

16 

(13.83,  14.64) 

107 

Total 

240 

129 

48 

417 

Table 

5.5. 

Population  tables 

considered  in 

this  simulation  study. 

e 

A = .5 

B 

A = .1 

-.576 

.499096  .343107  .157796  -.104 

.336749  .390702 

.273549 

0.000 

.310307  .379386  .310307  0. 

.304202  .391596 

.304202 

.576 

.157796  .343107  .499096  .104 

.273549  .390702 

.336749 

B 

A = .5 

B 

A = .1 

-.384 

.499096 

.343107 

.157796 

-.068 

.338892 

.386417 

.273549 

.192 

.309356 

.381289 

.309356 

.034 

.306354 

.387292 

.306354 

.192 

.309356 

.381289 

.309356 

.034 

.306354 

.387292 

.306354 

Note: 

tables 

The  two 
pertain 

tables  at 
to  case  2. 

the  top 

pertain  to 

case  1 . 

The  remaining 
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Table  5. 6. a.  Computed  statistics  for  case  1 (A  = .5) 


Sample 

Row 

Statistics 

size  N 

1 

2 

3 

Parameter 

Bi 

-.576 

0 

.576 

6 

30 

-.664 

.005 

.658 

i 

150 

-.598 

.0005 

.597 

★ 

30 

-.688 

.003 

.684 

6i 

150 

-.598 

.0002 

.598 

30 

.217 

.162 

.208 

ni 

150 

.027 

.025 

.030 

-v* 

30 

.195 

.103 

.185 

ni 

150 

.027 

.024 

.029 

30 

.128 

.113 

.128 

ni 

150 

.026 

.025 

.026 

A 

30 

.006 

.005 

.006 

ni 

150 

.0011 

.0011 

.0011 

A 

30 

.013 

.008 

.013 

n,- 

150 

.0013 

.0012 

.0014 

30 

.012 

.005 

.012 

a~* 

ni 

150 

.0012 

.0011 

.0014 

r'i/ni 

30 

.900 

.637 

.890 

150 

.987 

.958 

.980 

Note:  The  number  of  times  the  estimated  parameters  were  out  of  order  is 

257(172)  for  the  sample  size  N = 30  (N  = 150),  i.e.,  P(B-,^MBo)  = 

.743  (.828)  for  N = 30  (150).  1 2 3 
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Table  5 . 6 . b . 

Computed  statistics  for  case  1 

(A  = .1). 

Statistics 

Samp! e 
size  N 

1 

Row 

2 

3 

Parameter 

ei 

-.104 

0 

.104 

30 

-.107 

.003 

.104 

ei 

150 

-.111 

-.001 

.112 

★ 

30 

-.208 

-.010 

.218 

gi 

150 

-.140 

-.003 

.140 

30 

.154 

.144 

.151 

ni 

150 

.022 

.024 

.024 

30 

.099 

.031 

.088 

ni 

150 

.016 

.008 

.017 

30 

.110 

.110 

.110 

ni 

150 

.022 

.022 

.022 

30 

.005 

.005 

.005 

ni 

150 

.001 

.001 

.001 

A 

30 

.008 

.007 

.008 

ni 

150 

.010 

.001 

.001 

^ * 

30 

.007 

.002 

.006 

ni 

150 

.001 

.0004 

.001 

30 

.639 

.218 

.589 

Vi 

150 

.700 

.349 

.703 

Note:  The  number  of  times  the  estimated  parameters  were  out  of  order  is 

691  (595)  for  the  sample  size  N = 30  (N  = 150),  i.e.,  P(Bi<6o<§o) 

= .309  (.405)  for  N = 30  (150).  1_  2“  3 


147 


Table  5. 7. a. 

Computed  statistics  for  case  2 

(A  = .5). 

Statistics 

Sampl e 
size  N 

1 

Row 

2 

3 

Parameter 

-.384 

.193 

.193 

30 

-.439 

.226 

.213 

ei 

150 

-.403 

.197 

.206 

★ 

30 

-.461 

.118 

.343 

3i 

150 

-.402 

.146 

.256 

A 

30 

.191 

.154 

.160 

^i 

150 

.025 

.025 

.024 

30 

.165 

.061 

.109 

"l 

150 

.024 

.015 

.017 

30 

.122 

.111 

.111 

^i 

150 

.024 

.022 

.022 

Q/\ 

30 

.005 

.005 

.005 

ni 

150 

.001 

.001 

.001 

/A 

30 

.011 

.008 

.008 

ni 

150 

.001 

.0007 

.001 

o~* 

30 

.010 

.003 

.008 

150 

.001 

.0007 

.001 

/v*  /\ 

Vni 

30 

.863 

.397 

.680 

150 

.980 

.591 

.706 

Note:  The  number  of  times  the  estimated  parameters  were  out  of  order  is 

565  (473)  for  the  sample  size  N = 30  (N  = 150),  i.e.,  P(6,<MfU  = 

.435  (.527)  for  N = 30  (150).  1 2 3 
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Table  5.7.b.  Computed  statistics  for  case  2 (A  = .1). 


Sample 

Row 

Statistics 

size  N 

1 

2 

3 

Parameter 

*1 

-.068 

.034 

.034 

6 

30 

-.067 

.041 

.027 

i 

150 

-.069 

.034 

.034 

* 

30 

-.168 

-.001 

.169 

Pi 

150 

-.104 

.008 

.096 

30 

.152 

.140 

.148 

*1 

150 

.022 

.024 

.023 

30 

.095 

.027 

.082 

ni 

150 

.015 

.007 

.015 

30 

.109 

.109 

.109 

ni 

150 

.022 

.022 

.022 

30 

.005 

.005 

.005 

*1 

150 

.001 

.001 

.001 

A 

30 

.008 

.007 

.008 

*1 

150 

.001 

.001 

.001 

30 

.007 

.001 

.006 

ni 

150 

.001 

.0003 

.001 

. ★ - 

r'i/T'i 

30 

.624 

.194 

.555 

150 

.696 

.289 

.640 

Note:  The  number  of  times  the  estimated  parameters  were  out  of  order  is 

723  (712)  for  the  sample  size  N = 30  (150),  i.e.,  P(0,<|L<(L)  = 

. 273 ( . 288 ) for  N = 30  (150).  1 1 6 
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Table  5.8.  Relative  MSE  for  different  shifts,  sample  sizes,  and 
category  distributions. 


Shift 

Sample  size 

The  three  categories  The  last  two  categories 

have  different  distributions  are  identical 

A = .5  A = .1  A = .5  A = .1 

N = 30 

.824  .487  .663  .464 

N = 150 


.976 


.580 


.760 


.537 


CHAPTER  SIX 
SUMMARY 

In  this  dissertation,  we  discussed  association  models  for  cross- 
classifications having  ordered  categories.  If  the  categorical 
variables  are  treated  as  if  they  were  nominal  in  scale,  Bishop  et  al . 
(1975)  give  a comprehensive  treatment.  However,  if  some  or  all  of  the 
variables  have  a natural  order,  then  it  is  possible  to  construct 
simpler  models  than  if  the  variables  are  treated  as  if  they  were 
nominal  in  scale.  In  the  models  designed  for  ordinal  data,  one  can 
either  assign  or  estimate  scores  of  the  levels  of  an  ordinal  variable. 

In  Chapter  Two,  we  reviewed  the  models  for  two-way  tables  intro- 
duced by  Goodman  (1979).  However,  we  presented  them  in  a different 
form  to  clarify  how  they  relate  to  models  designed  for  nominal 
variables.  In  Chapter  Three,  we  extended  them  to  multi-way  cross- 
classifications. We  discussed  their  interpretation  using  examples  and 
we  also  gave  a list  of  available  procedures  for  fitting  them.  Unlike 
in  models  designed  for  nominal  variables,  the  user  has  a rich  variety 
of  models  to  choose  from  if  the  natural  ordering  of  the  variables  is 
taken  into  account.  However,  this  great  variety  of  models  can  make 
model  selection  more  difficult.  Hence  future  work  should  be  directed 
at  finding  strategies  in  selecting  models. 

In  Chapter  Four,  we  proposed  statistics  for  testing  independence 
in  a two-way  table.  Through  a simulation  study,  it  appeared  that  the 
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new  statistics  are  more  powerful  for  detecting  certain  associations 
between  the  variables  than  the  conventional  chi-square  statistics  for 
testing  independence,  whenever  there  is  a trend  in  the  data.  Haberman 
(1981)  proposes  different  statistics  for  testing  independence.  However, 
we  did  not  compare  his  statistics  to  ours.  Thus,  it  will  be  desirable 
in  the  future  to  conduct  a simulation  study  to  compare  all  the 
statistics  designed  to  test  for  independence  between  the  rows  and 
columns  of  a two-way  table. 

Whenever  scores  cannot  be  assigned  to  the  levels  of  an  ordinal 
variable,  it  is  possible  to  estimate  them.  However,  the  estimated 
scores  need  not  be  monotonic.  To  overcome  this  shortcoming,  in 
Chapter  Five  we  proposed  methods  that  produce  isotonic  estimates  of 
the  association  parameters  in  the  R model  (2.2)  and  the  RC  model  (2.5). 

A simulation  study  showed  that  if  the  association  parameters  are  truly 
ordered  in  the  R model,  the  proposed  estimates  perform  better  than  the 
regular  estimates  in  mean  square  error.  For  the  case  of  the  R model 
and  RC  model,  we  conjecture  that  the  proposed  estimates  are  ML  estimates, 
so  that  future  work  should  be  conducted  in  that  direction.  It  will  be 
also  of  interest,  in  the  future,  to  extend  the  methods  above  to  multi- 
way cross-classifications. 


APPENDIX  1 

ASYMPTOTIC  NORMALITY  OF  THE  ML  ESTIMATES  OF  MODEL  PARAMETERS 

Let  us  consider  the  cell  proportions  {p..}  as  functions  of 

1 J 

6_  = (0-|,...»0  ),  a vector  valued  parameter  belonging  to  a set  0. 

The  true  value  e®  is  supposed  to  be  an  interior  point  of  0.  We 
therefore  make  the  following  assumptions: 

(i)  Pjj(0)  * Pij(j>‘)  for  at  least  one  (i,j)  when  .0  * _e'  (weak 
identifiabil ity  condition). 

3 P,-  i 

(ii)  All  the  partial  derivatives  {■---■ are  continuous  at  the  true 

30 

0 s 
value  e . 

(iii)  The  information  matrix  1(0)  is  nonsingular  at  0 =6°. 

Then  Q_  has  an  asymptotic  normal  distribution  (Rao.  C.R.,  1972,  p.  295). 

To  prove  the  asymptotic  normality  of  the  ML  estimates  of  the  model 
parameters,  we  have  to  verify  the  above  assumptions.  We  will  verify 
them  for  the  1 inear-by-1 inear  association  model.  It  can  be  done  in  a 
similar  way  for  the  other  association  models  for  two-way  tables.  The 
1 inear-by-1 inear  association  model  is: 

P-jj(9)  = exp(y+X^+Xj+3(ul.-u)(v.-v’))  l£i£r,  l<j<c,  where 

^ Y v v y 

ZXi  = EXj  = 0.  In  this  case  the  vector-parameter  0^  = (y,Xp...,X*  ^ ,x| 

» • • • » 6 ) 

(i)  Assume  p . .(e)  = p,  .(0')  l<iSr,  l<j<c.  Hence  log  p.  .(e)  = log  p. . (o' ) 

' J IJ  I J "I  J 
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>-y+A*+AY+  B(u.|.-u)(Vj-v)  = y + A *+AjY+  3'  (u^ -u)  (v^-v) . 


(A.  1) 


I 

If  we  sum  over  all  (i,j)  in  (A.l)  we  obtain  that  y = y . Now  sunming 

X 1 X 

over  j,  we  obtain  A.  = A^  l^i^r,  and  summing  over  i we  obtain 
Y ' Y 

Xj  - Aj  lij^c.  Therefore,  both  vector  parameters  £ and  _0‘  are  equal. 


9pi  • 

(11)  nr1  = pij-  3Xx 


3 Pi  j 3Pii  9 p.  . 

1 = 6 .p. .,  = 6 .p. . and  — ^ = 


'ri-ij’  3XV 


sru 
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Z(u1-u)(vj-v)pij. 

where  6kk,  is  the  Kronecker  delta.  Note  that  all  the  partial  derivatives 
are  continuous  at  the  true  value  0°.  Note  also  that  the  partial 
derivatives  of  higher  order  will  be  continuous  at  e_  = e°.  In  addition  to 
that,  they  are  bounded  around  the  true  value  6°. 

(iii)  We  will  assume  that  the  information  matrix  is  nonsingular  at  0°. 


APPENDIX  2 


The  listing  of  the  SAS  computer  program  used  to  fit  the  constrained 
R Model . 
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**************  ******************  ******* 


* * 

* THIS  PROGRAM  FITS  THE  CONSTRAINED  * 

* E-MODEL  TO  AN  RXC  CROSS-CLASSI-  * 

* FICATION  USING  THE  METHOD  OF  * 

* SECTION  2 OF  CHAPTER  5.  * 

* * 


*************************************** 

DATA  FREQ; 

INPUT  Z1-Z4; 

CARDS; 

INPUT  THE  DATA  IN  A TABLE  FORM 
PEOC  MATRIX; 

WE  INITIALIZE  THE  PARAMETERS 


FETCH  FF  DAT A= FREQ; 

F=FF ' ; 

PRINT  F; 

fi=NSO  W (F)  ;C=NCOL (F)  ; 

NR=S+C; NC=R+C- 1 ; 

NP=SHAPE (F,NB) ; N=NP* ; 

NI=N; 

NN=DI AG  ( N)  ;L=LOG  ( II)  ;LI=L; 

XR= J (NE#E-1f0)  ;U=J(E,1)  ;V=J(C,1)  ; 

UT-J  (R,  1 ,0) ; UR=J  (H- 1,1,0)  ;ORB=J(R-1,  1,0)  ; 

WE  DETERMINE  THE  DESIGN  MATRIX 


X=J  (NR,  NC,  0)  ; 

DO  1=1  TO  NIJ;X  (I,  1)  =1  ;END; 

DO  J=1  TO  E- 1 ; DO  1=1  TO  C; 

X(I+(J-1)  *C,J*1)  =1;  END;  END; 

DO  J=1  TO  R- 1 ; DO  1=1  TO  C; 

X(I+(E-1)  +C,J+1)  =-1  ;£ND;END; 

DO  J=1  TO  C-1  ; DO  1=1  TO  S; 

X ( (I-  1)  ♦ C + J,  J + R)  =1  ; END; END; 

DO  J=1  TO  C-1;  DO  1=1  TO  S; 

X (I*C, J+R) =-1 ; END; END; 

DO  J=1  TO  C;V  (J,  1)  =J- (C+1)  #/2;END; 

MAIN  PROGRAM 

LINK  ROKEF; 

If  DD=0  THEN  GO  TO  S3; 

ELSE  LINK  CEDRjLINK  ROKEF; 

S3:  RETURN; 

THIS  PROCEDURE  FITS  THE  R MODEL 
TO  THE  RXC  TABLE 
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ROWEF: 

ELSE  DD=1 ; 

RETURN; 

ORDR: 

RETURN; 

DO  J=1  TC  R-1  ; DO  1=1  TO  C; 

XE  (1+  (J-1)  *C,  J)=V  (I#  1)  ; END;  END; 

DO  J=1  TO  fi-1  ; DO  1=1  TC  C; 

XE  (1+  (R-1)  *C.  J)  =-V  (1,1)  ; END;  END; 

XI  =X | | XR; 

NP=SH  APE  ( F,  NR)  ;N=NP'  ;L=LOG  (N)  ; 

NN=DI AS  (N) ; H1  = X1 ' *NN*X1 ; 

R1  = X1  ' * N N *L  ; B 1 = (INV  (HI)  ) *Rl ; 

M 1=EXP (X1+B1) ; 

DO  K=  1 TO  5; 

MM1=DI AG  (Ml)  ; Q1  = Xl'  * (N-M1)  ;H1=X1 • *HN1*X1; 

BE 1=B 1+ (INV(HI) ) *Q1 ; 

M1  = EXP(X1*BB1)  ; B1  = BB  1 ; END; 

GSQ=2#NI • * (LI- LOS (Ml) ) ; 

XSQ=  ( Nl-iil)  • * I NV  ( M K 1 ) * (NI-Ml)  ; 

PRINT  Ml  GSQ  XSQ  ; 

DO  1=1  TC  R-1  ; URE (I,  1)  =B1  (I+R  *C-  1 , 1)  ; END  ; 

SM  U=- SUM  (URR)  ; UT=URR//SMU; 

U=UT ; PEI  NT  U; 

UH(  1;R-1r1)=UT(1;  R-  1 # 1 ) — UT  (2  ; R # 1 ) ; 

IF  MAX(Ufi)  < .00  1 THEN  DD=0; 

ELSE  DD=1 ; 

RETURN; 

THIS  PROCEDURE  TRANSFORMS  THE  TABLE 
ACCORDING  TO  EQUATION  (5.5) 

ORDR: 

FR=F(  ,+) ; K=1 ; 

SSI:  KK=K  + 1;XX=0.  0;  XX  = U (KK,  1) -U(K,  1)  ; 

IF  XX  > .001  THEN  GO  TO  SS2 ; 

DO  J=1  TO  C;Y=(F(K,J)  +F(KK, J)  ) #/  (FE  (K,  1)  +FR  (KK,  1)  ) ; 
F (K,  J ) = Y*FR (K,  1 ) ;F(KK,J)  =Y*FR(KK,1) ; EN  D ; 

SS2:  IF  K=R— 1 THEN  GO  TO  SS3;K=K+1;GO  TO  SSI; 

SS3:  PRINT  F ; RETURN ; 
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